Matheuristic approaches to the green sequencing and routing problem

Abstract

The paper addresses the green sequencing and routing problem, which consists in determining the best sequence of locations to visit within a warehouse for storing and/or retrieval operations, using a fleet composed of both electric vehicles, e.g., equipped with a lithium-ion battery, and conventional vehicles, i.e., with an internal combustion engine. We present a Mixed-Integer Linear Programming formulation to the problem and propose two matheuristics based on suitable decompositions of the mathematical formulation. The two matheuristics have been tested on a pool of small-medium size instances and their performance has been compared to the one of a third matheuristic, previously proposed for the case of conventional vehicles only and here suitable extended to deal with the green aspects of the problem. The performed analysis allowed one to identify the most promising matheuristic in terms of some standard computational indicators, i.e., computing time and percentage optimality gap, as well as in terms of some qualitative aspects of the solutions agreed with a reference company. Such a most promising algorithm has then been further tested to gather some technical insights on what makes the problem hard to solve, as well as to outline some managerial insights. Moreover, its performance has been tested on a pool of real instances comprising ordinary days (with a usual amount of operations to perform) and extremely busy days, showing its efficacy and efficiency also in the considered real application context.

1 Introduction

Warehouses are a basic component of any supply chain. In many of them, especially in Western Europe, workers walk or drive through the warehouse to perform picking and/or put-away operations (van Gils et al. 2018). The former concern the movement of items from the storage locations towards the output point of the warehouse to meet customer demands, while the latter concern the movement of items from the input point of the warehouse towards the storage area, to store the items in the assigned storage locations. Picking and put-away are recognized as the most labour and time-consuming operations within the warehouse, and their efficient planning plays a relevant role in both improving productivity and reducing the overall operational cost.

The problem addressing the issue of determining the most efficient sequence of operations to move items within the warehouse, by performing order picking and/or put-away operations, is known in the literature as Sequencing and Routing Problem (SRP). Typically some additional constraints, usually related to the considered application context, have to be satisfied as well, while minimizing the material handling cost or some travel measures, e.g., expressed in time or distance travelled by the workers (van Gils et al. 2018).

Considerations about environmental degradation, greenhouse gas (GHG) emissions, fossil fuel depletion and global warming have led many companies to increase their attention towards objectives such as sustainability and environment-friendly issues in warehouse management, in addition to the traditional operational and economic ones (Bartolini et al. 2019).

The Green SRP (GSRP) is thus emerging as a new topic of research. It is a generalization of SRP where also electric vehicles (EVs) perform picking and put-away operations. The use of EVs reduces GHG emissions promoting sustainability in logistics, reduces long-term management costs and improves healthiness for workers (e.g., reduced noise and better local air quality) as well as the company’s ability to quickly and flexibly respond to market changes (Jiao et al. 2021a, b; Bartolini et al. 2019). Additionally, an environment-friendly approach increases company reputation and brand fidelity, since customer awareness has changed towards improved sustainability. Nevertheless, considering EVs increases the problem complexity since specific activities, like the scheduling of recharging periods and the limited autonomy of such vehicles, have to be considered when planning picking and put-away operations.

The number of contributions on GSRP available in the literature is still quite limited, highlighting the novelty of the topic (Carli et al. 2020a). Also, very few case studies and realistic applications can be found on the subject, as outlined in Bartolini et al. (2019). There are, however, interesting contributions discussing battery management issues in GSRP. Some of them consider EVs performing operations as long as their battery charge is available, being then unavailable till the full recharge of the batteries (Ene et al. 2016; Khoei et al. 2022). Others, instead, apply the so-called battery swapping, i.e., the replacement of a low battery with a fully charged one when needed (Basso et al. 2019a; Carli et al. 2020a, b). A more detailed overview of this subject will be provided in the following section.

Recently, Lanza et al. (2022d) addressed an SRP where conventional vehicles, i.e., with an internal combustion engine, are used to perform picking and put-away operations in a given time horizon. The work has been motivated by a real application involving a large logistics site of an Italian company located in Tuscany. The involved warehouse is larger than 10,000 $\hbox {m}^2$, has a rectangular internal layout composed of narrow storage aisles and wide cross aisles, and comprises almost entirely storage areas. Thus, the distance travelled to perform operations is very large. More than 300 different types of tissue products for domestic and sanitary use are produced in this site, which is characterized by all-day long operations (i.e., 24 h per day) and by a high product rotation index (i.e., more than 1,000 product units are moved per day). After that work, a preliminary study has been carried out by the same authors to investigate the green extension of the above-mentioned SRP, where some of the vehicles are electric and equipped with a lithium-ion battery (Lanza et al. 2022c). This kind of battery is considered the most promising for the near future by the majority of the literary sources, for its high energy efficiency, long lifespan and the possibility to be frequently charged, even partially, during short break times between operations, without being removed (Andwari et al. 2017). Additionally, such a battery may be slightly recharged when the vehicle brakes, if the vehicle is equipped with the regenerative braking recovery mechanism that converts the kinetic energy produced into electrical energy. Indeed, this is also the technology adopted in the warehouse of the Italian company mentioned before. Specifically, in Lanza et al. (2022c) the authors proposed an extension of the mathematical model in Lanza et al. (2022d), aiming at planning picking and put-away operations for both conventional and EVs, and, for the latter, also scheduling break times to recharge the battery at the available charging station, by exploiting the partial recharging feature proper of lithium-ion batteries. This was the first modelling attempt to apply the partial recharging policy while scheduling both picking and put-away operations in a labour-intensive warehouse, instead of considering full battery recharge or battery swapping policies as previously modelled in the literature. The authors also generalized the matheuristic approach to SRP presented in Lanza et al. (2022d), which is based on a time horizon decomposition, to deal with the green features of the extended problem.

In this paper, we further enhance the study on GSRP by extending the results in Lanza et al. (2022c). Firstly, we generalize the mathematical model in Lanza et al. (2022c) regarding the green aspects. Specifically, we extend the objective function by adding a further soft objective aimed at avoiding unnecessary periods of charge of the EVs in the warehouse, i.e., those periods where an EV is charged even though it does not have any forthcoming planned operation to perform. Unnecessary recharge periods may raise the cost of warehouse management, uselessly consume resources within the warehouse, and increase the environmental impact (Bartolini et al. 2019). Moreover, we add to the model constraints of practical impact regarding the management of EVs. Precisely, we introduce battery management constraints which prescribe a minimum state of charge of the battery of EVs at the end of the considered planning horizon. Such constraints are particularly relevant for those warehouses performing all-day-long operations, as happens in the case study addressed, where it is required that, at the beginning of the next time horizon, the state of charge of each EV be sufficient to perform some basic operations, rather than being completely discharged. Also, we introduce capacity constraints for the charging station present in the warehouse, which impose a limit on the number of EVs that may be charged simultaneously. All such additional features greatly increase the difficulty of GSRP for the version studied in Lanza et al. (2022c). Notice that both SRP and GSRP are NP-hard, being generalizations of the Capacitated Vehicle Routing Problem (Masae et al. 2020), and that the additional green aspects addressed in this work render the problem much more adherent to the reality of all-day working warehouses, but even harder to solve.

A consequence is that real-size instances of GSRP, such as the ones related to the considered case study, can not be directly addressed by solving the proposed mathematical formulation via a state-of-the-art commercial solver. For this motivation, starting from the presented mathematical formulation, we propose and study three matheuristic approaches, to compute solutions of good quality in a reasonable amount of time. Precisely, firstly we extend the time horizon decomposition matheuristic already proposed in Lanza et al. (2022c), to deal with the new green features of the problem previously described. Then, we propose two new matheuristic algorithms to tackle GSRP. The first algorithm performs an activity decomposition and it consists of two phases: in the first phase put-away activities are scheduled, while in the second phase picking activities are planned given the partial schedule obtained in the first phase. The second algorithm, instead, combines the features of the above-mentioned activity decomposition matheuristic with the time horizon decomposition approach, by performing both a time and an activity decomposition.

We have firstly analysed the computational performance of the three algorithms mentioned above on a set of instances of small-medium size based on real data (a larger set compared to the test bed used in Lanza et al. 2022c), to provide a comparative evaluation of the three matheuristics in terms of computational performance, i.e., computing time and optimality percentage gap, as well as in terms of some qualitative indicators suggested by the management of the reference warehouse (i.e., the one considered in the reported case study). The evaluation has shown that the algorithm combining activity and time decomposition is particularly efficient in the tested instances and produces solutions of good quality in terms of the considered indicators. Such an algorithm has therefore been used to gather both some technical insights on which aspects of GSRP make the resolution process particularly complex, and also some managerial insights focusing on green aspects such as the number of used EVs and the number of available charging stations. Finally, the algorithm combining activity and time decomposition has been used to solve a set of real instances shared by the reference company, which comprises 5 ordinary and 5 extremely busy working days, showing its efficiency and efficacy also in the real application context addressed.

The paper is organized as follows. Section 2 is devoted to the literature and reviews the main results in the area of SRP (in Sect. 2.1) and the relevant yet sparse literature on GSRP (in Sect. 2.2). Section 3 describes the GSRP addressed in this paper in a detailed way. The mathematical model is presented in Sect. 4, by firstly introducing the main notation used to formulate GSRP (Sect. 4.1), then discussing battery-related issues and the energy consumption model (Sect. 4.2), and finally presenting the objective function and those mathematical constraints which are peculiar to GSRP (Sect. 4.3). The overall mathematical formulation of GSRP can be obtained by adding the set of constraints reported in Appendix A, which model the aspects of the problem with only conventional vehicles. Section 5 is devoted to the matheuristic approaches. The time horizon decomposition matheuristic, tailored to GSRP, is presented in Sect. 5.1. Then, the two new algorithms are described: the one based on activity decomposition (Sect. 5.2) and the one combining time horizon and activity decomposition (Sect. 5.3). The results of the wide computational experimentation are presented in Sect. 6: the performance of the proposed matheuristics on a set of instances of small-medium size is evaluated in Sect. 6.4.1, some technical and managerial insights are provided in Sect. 6.4.2, while the most promising algorithmic approach, i.e., the one combining time horizon and activity decomposition, is tested on a set of real instances in Sect. 6.5. Finally, Sect. 7 concludes the paper and identifies some future research directions, which focus on different battery management policies, such as the battery swapping, as well as the integration of the GSRP with another crucial operational problem typical of warehouse management, i.e., the storage location assignment problem.

2 Literature review

We summarize here some of the most recent and influential contributions available in the literature related to SRP (in Sect. 2.1) and GSRP (in Sect. 2.2), then discuss the positioning of our contribution relative to the literature (in Sect. 2.3).

2.1 SRP

The literature dealing with picker-to-part warehouse systems, i.e., systems where workers walk or drive through the warehouse to pick up and/or put away items, has focused almost exclusively on designing picking routes, making contributions focused on the combination of both picking and put-away much more scarce.

The most recent contributions for picking activities only focus on realistic aspects such as particular layout designs (Mowrey and Parikh 2014; Scholz et al. 2016; Boysen et al. 2017; Weidinger et al. 2019; Briant et al. 2020), congestion issues (Pan and Wu 2012; Chen et al. 2013, 2016), workers comfort (Grosse et al. 2015) and dynamic modification of list of operations (Lu et al. 2016; De Santis et al. 2018). As opposed, Gómez-Montoya et al. (2020) is the only contribution addressing exclusively a put-away SRP. Integrated schedule planning of storage and picking operations certainly define difficult larger problems to tackle, providing however the opportunity to assign resources more efficiently and have better performances from a practical point of view. Contributions in which joint storage and retrieval operations are considered are Wruck et al. (2013), in which SRP is modelled as a multi-objective minimization model for the single-worker case, in both static and dynamic settings; Schrotenboer et al. (2017), in which SRP is modelled as a variant of the travelling salesman problem for both the single-worker case and the multiple-worker case; and Ballestín et al. (2013), in which SRP is modelled as a project scheduling problem, in both static and dynamic settings. A crucial aspect that may influence the system’s performance is the layout of the warehouse. Pohl et al. (2009a) consider the three most common rectangular layouts, while Pohl et al. (2009b) and Gue et al. (2012) address uncommon layouts such as the Flying-V and Inverted-V aisles design. For additional contributions on SRP we refer to relevant surveys such as van Gils et al. (2018) and Masae et al. (2020).

2.2 GSRP

Sustainability, as a general theme, has received increased consideration in supply chains in the last decades. The green supply chain management literature covers various topics and deals with miscellaneous aspects, including green logistics, reverse logistics, green building, product life cycle assessment, or environmental and operational performance (de Oliveira et al. 2018; Waltho et al. 2019; Bartolini et al. 2019). The term green warehousing has been introduced to specifically denote the integration and implementation of environment-friendly operations in a warehouse, to minimize energy consumption, energy cost and GHG emissions (Bartolini et al. 2019). Among the different macro-themes in green warehousing, energy saving in warehousing focuses on the optimization of operations performed by both fixed (e.g., conveyor systems, automated cranes) and mobile material handling equipment (e.g., electric forklift, electric automated vehicles) from an energy-efficiency perspective, thus considering their limited autonomy and the need of recharging periods for their batteries (Bartolini et al. 2019). This is the case of GSRP.

In the following, we shall overview some contributions to energy saving in warehousing by considering manual picker-to-part systems, i.e., systems where operations are performed by workers walking or driving along the aisles of the warehouse. The reader interested in a broader discussion on green supply chain management is referred to the surveys of Srivastava (2007); Biel and Glock (2016); de Oliveira et al. (2018); Bartolini et al. (2019); Waltho et al. (2019); Bänsch et al. (2021).

A recent but well-studied line of research analyzes battery issues for EVs (Thangavel et al. 2023; Habib et al. 2023). Private and public transportation is mostly considered in the literature (see the review of Andwari et al. 2017). Nevertheless, some contributions in warehousing are also available. Without making routing decisions or planning battery charging periods, but rather considering average estimations on the energy consumption, the focus of these studies is to enable managers to evaluate optimal infrastructure requirements, e.g., forklift type, type of battery, storage configuration, warehouse layouts or picking policy (see for instance Boenzi et al. 2015; Facchini et al. 2016, 2018; Lewczuk et al. 2021; Jiao et al. 2021b, a).

Focusing on routing aspects, instead, the target of most proposed models is to improve the environmental performance of the system by minimizing the energy consumption of the vehicles while performing operations. The estimation of the energy consumption may be a function of the time travelled by vehicles, therefore considering so-called horizontal moves only (see for instance Vivaldini et al. 2010; Confessore et al. 2013; Anand et al. 2014; Rojanapitoon and Teeravaraprug 2018; Burinskiene et al. 2018), or it may also include the energy consumed by a vehicle to lift items, while the energy to lower items is usually assumed negligible instead. Ene et al. (2016) propose a genetic algorithm to solve an order batching and routing problem for picking operations only, in which the energy consumption evaluation is based on the energy consumed per time unit, without considering load, acceleration, deceleration, and friction forces. Khoei et al. (2022) consider electric forklifts equipped with an elevating platform from which standing operators perform picking operations from racks. No due time for bringing the ordered items to the depot location is considered. The horizontal energy consumption is calculated by considering friction forces, acceleration and deceleration, while the vertical energy consumption is calculated based on the height the platform needs to achieve to let the operator pick an item, and on the number of items already picked by the operator and collected on it. A mixed-integer programming formulation is proposed as well as dynamic programming and heuristic approaches.

The two above-mentioned contributions do not explicitly consider battery management in their formulation, rather it is sometimes only assumed that an EV operates until the state of charge of its battery is higher than a threshold. Then, it becomes unavailable till the full recharge of the battery. Other contributions, instead, consider different battery management such as the replacement of a low battery with a fully charged one when needed, the so-called battery swapping. If, on the one hand, this policy allows it to have EVs always be available to operate, on the other hand, it requires an efficient schedule to recharge the batteries, to ensure a significant number of exploitable batteries for swapping (Andwari et al. 2017). Lee et al. (2022) propose an energy-aware dynamic model that is capable of establishing picking/put-away routing for a fleet of electric forklifts simultaneously with battery replacement operations at the charging area (battery charging schedules are not addressed). The authors propose a dynamic algorithm whose performance is evaluated by using real energy and vehicle operational parameters. Basso et al. (2019a) focus instead on the complementary problem of developing a schedule for charging batteries in a battery centre of a warehouse, where electric stackers are used to move goods. The problem involves assigning low-charge batteries to specific chargers and determining charging start times, considering both periods where inexpensive but limited photovoltaic solar energy is available and periods of varying prices of grid-available electric energy. The authors propose a scheduling jobs-in-parallel-machines formulation to minimize both the energy cost and the delay to fulfil the battery. The operations of the EVs outside the battery centre are not addressed. Carli et al. (2020a) develop an integer programming formulation to jointly determine the activities of a fleet of electrically powered forklifts and the charging schedules for batteries in a labour-intensive warehouse. During their activities, the forklifts may replace empty batteries with fully charged ones, which are then recharged at charging stations. The energy consumption is assumed to be proportional to the activity duration, while the energy costs are time-dependent. Precisely, two tariffs are considered, i.e., on-peak and off-peak. The goal of the model is to minimize the sum of the penalty cost related to the makespan over all the activities, and the total electricity cost for charging batteries. In Carli et al. (2020b), the authors introduce priorities among activities in the problem formulation and define a two-step optimization approach: the first step defines a high-level schedule of activities to be executed by forklifts following the given job priorities, but neglecting battery replacements; the second step, instead, determines the detailed schedules of the activities and the optimal battery charging strategies.

2.3 Positioning our contribution with respect to the literature

The reported overview highlights the two major battery management policies which have been addressed in the literature in the context of GSRP. That is to say, EVs perform operations as long as battery charge is available, being then either unavailable for operations for a long time until the full recharge of their batteries (e.g., an overnight recharge), or just a few minutes for the battery swap. Notice that, in the latter case, planning pick-up and put-away operations at the same time as battery recharge operations has not yet been taken into account.

In this study, firstly we consider a different policy related to battery recharge, consisting of the possibility to frequently recharge the lithium-ion batteries equipped on the EVs, without removing them from the vehicles, even during short break times between pick-up and put-away operations. Such a feature was investigated preliminarily in Lanza et al. (2022c), but here it is enhanced by considering additional energy constraints, related to a minimum level of charge of the EV battery at the end of the planning horizon, and to a maximum number of vehicles that may simultaneously be recharged, and also considering an additional term in the objective function, related to unnecessary periods of charge of the EVs when no forthcoming operation is planned. Such enhancements render the studied problem still more adherent to the reality of warehouses such as the one analysed in the considered case study, but also much more complex. This has motivated the design of clever matheuristics, based on the proposed mathematical formulation via decomposition, to efficiently solve also instances of large size.

We also point out that, even though lithium-ion batteries are recognized as a promising technology by the majority of literary sources, thanks to their high energy efficiency and long-lasting lifespan, the only few contributions available in the literature focus on a correct estimation of their energy usage, or on a comparison with more dated methods (see for instance Jiao et al. 2021b, a), rather than addressing the joint planning of operations within the warehouses and recharge periods, as we do.

Indeed, lithium-ion battery technology and related partial recharging have been considered in the larger area of green logistics and last-mile deliveries (for a review of green vehicle routing problems see Lin et al. 2014; Moghdani et al. 2021; Asghari and Al-e 2021). These kinds of problems, however, are different in that they do not consider demand rates, multiple vehicles per route, and multiple visits to the same location, as workers within a warehouse normally do.

The attributes of the aforementioned papers in Sect. 2.2 concerning routing problems of EVs in warehouses have been succinctly summarized in Table 1, which also outlines the distinctive features of the problem addressed in this study.

Table 1 Recent studies on EV routing problems

3 The problem

The kind of warehouse we consider is composed of two disjoint areas. The first area is a transit zone connecting the input points, where items wait to be stored, to the storage area. In the case study addressed, for example, the input points are end-of-line conveyor belts. The second area, instead, is devoted to the stock. This storage area is organized in separate departments where items are stocked in storage locations, homogeneously for the type of product. For example, in the considered case study the storage locations are stacks, where tissue products for sanitary and domestic use are stocked homogeneously. The output point of the warehouse, called the collection area, is located in this area. In the collection area, retrieved items are gathered to establish order integrity before loading trucks and shipping them. Storage locations may have different capacities, and both the input points and the collection area are capacitated.

During a given time horizon, items of different product types are placed on the input points and require transportation to their preassigned storage locations. These items are called incoming items. At the same time, items of different types need to be picked from their storage locations and transported to the collection area, to meet customer demands. Such items are called outgoing items. Each incoming item is available at some input point at a known availability date, while each outgoing item must reach the collection area before a known due date. Types and numbers of both incoming and outgoing items are known in advance and provided in a storing list and a shipping list, respectively.

Picking and put-away of the items are performed by capacitated vehicles belonging to two fleets having different characteristics, which will be referred to as F1 and F2 in the following. In the considered case study, F1 are autonomous LGV vehicles, while F2 are human-driven forklifts. The routing of these two fleets of vehicles is restricted to only one of the above-described disjoint areas of the warehouse. Specifically, vehicles of type F1 can only move in the transit zone, whereas vehicles of type F2 can only move in the storage area. Vehicles of different fleets may exchange freights at specific capacitated points called collectors, which are located along the boundary of the two disjoint areas. Specifically, incoming items are picked up from the input points by a vehicle of type F1 and transported to one of the available collectors, where they are unloaded. From there, items are loaded by a vehicle of type F2, that moves them from the collector towards their preassigned storage locations. On the other hand, outgoing items are loaded from their storage locations by vehicles of type F2 and moved straightforwardly to the collection area.

A possible structure of the warehouse is shown in Fig. 1, where a rectangular layout organized in four departments is depicted. The figure reports the position of the input points (denoted with IP), of the collectors (denoted with C), and of the unique charging station (located in department 2). Different colours, i.e., light grey and dark grey, denote the areas of the warehouse where vehicles of type F1 and F2 can move, respectively. Additionally, possible itineraries for incoming and outgoing items are outlined with white arrows. Precisely, the figure shows, as an example, the itinerary of an incoming item released on one of the three depicted input points and transported to a storage location in Department 1, and the itinerary of an outgoing item picked from Department 4 and moved toward the collection area.

Fig. 1

A possible warehouse structure

In the green context under consideration, some vehicles of type F2 are electric. Specifically, inspired by the specific type of EVs used in the warehouse of the case study addressed, EVs equipped with a lithium-ion battery and with a regenerative braking recovery mechanism are considered. The battery is discharged when the vehicle moves or when it lifts items from the ground, while it is slightly recharged when the vehicle brakes, due to the regenerative braking mechanism mentioned earlier. Recharge operations are performed at a specific charging station, even during short break times between operations, without removing the battery from the vehicle. In recharging EVs, a specific requirement is that the level of charge of their battery be maintained within a given range, to ensure a long lifetime to the battery itself. Moreover, it is required that the fleet of EVs be not completely discharged at the end of the considered time horizon, but rather the amount of charge be adequate to be operative in the next time horizon. Application contexts where such a requirement is critical are warehouses that operate all day long (i.e., 24 h per day), as in the case study addressed. Requirements related to battery management are fully detailed in Sect. 4.2.

The goal of the studied problem is to plan the routing of all the vehicles to move the incoming items from the input points toward the preassigned storage locations, and the outgoing items from storage locations toward the collection area, with the primary aim of minimizing the total vehicle travel time within the warehouse. In achieving this, battery management and recharging operations must be planned for the EVs composing the fleet.

Indeed, some additional features have been considered during the optimization process. They involve the anticipation movement of some outgoing items towards the collection area concerning their planned due dates, given the high number of operations expected during the planning horizon. Also, a strict management policy for both picking and put-away operations is imposed, separately per product type. That is, storage locations have to be emptied/filled up one at a time following a prespecified order of precedence, due to the perishability of the products in the warehouse, like in the considered application context, with the consequent need to retrieve and ship first the items with the highest time of permanence within the warehouse. Finally, safety requirements for workers prescribe that no two vehicles may travel from the same location toward another same location at the same time. For more details on such complicating management constraints, already present in the version of the problem with only conventional vehicles, we refer to Lanza et al. (2022d).

We conclude by observing that both the SRP and the GSRP are generalizations of the Capacitated Vehicle Routing Problem, and therefore they are both NP-hard. The additional green aspects presented before, and described in more detail in the following section, indeed make the problem even harder to solve.

4 Modeling GSRP

In this section, firstly we introduce the main notation used to formulate GSRP (Sect. 4.1), disregarding however the specific notation related to the battery management, which is presented in Sect. 4.2. Section 4.2 also describes the energy consumption model used to estimate the energy consumed by the EVs within the warehouse. Finally, in Sect. 4.3 we state those constraints which regulate the energy behaviour of the EVs, since they are peculiar to GSRP, by referring to Appendix A for the additional constraints that complete the mathematical formulation of the problem (such more standard constraints, which are common to the version of the problem with only conventional vehicles, have already been presented in Lanza et al. 2022d).

4.1 Model notation

Let ${\mathcal {K}}$ be the set of product types to move in the given time horizon. We distinguish the subset of product types of the incoming items, ${\mathcal {K}}_{in}$, and the subset of product types of the outgoing items, ${\mathcal {K}}_{out}$. Let ${\mathcal {V}}^{1}$ and ${\mathcal {V}}^{2}$ be the sets of vehicles of type F1 and F2, respectively. Moreover, let ${\mathcal {V}}^{E} \subseteq {\mathcal {V}}^{2}$ denote the subset of the EVs.

We denote by ${\mathcal {G}}^P=({\mathcal {N}}^P, {\mathcal {A}}^P)$ the directed graph representing the physical network along which vehicles move. The set of nodes ${\mathcal {N}}^P$ comprises:

• the set ${\mathcal {S}}^k_{in}$ of the storage locations preassigned to the product types in ${\mathcal {K}}_{in}$, and the set ${\mathcal {S}}^k_{out}$ of the storage locations occupied by items of product types in ${\mathcal {K}}_{out}$ at the beginning of the time horizon;

• the parking areas for vehicles of type F1 and F2, denoted by $\omega ^{1}$ and $\omega ^{2}$, respectively;

• the set ${\mathcal {R}}$ of the input points;

• the set ${\mathcal {B}}$ of the collectors;

• the output point, or collection area, $\pi$;

• the charging station c.

The set ${\mathcal {A}}^P$ represents direct connections between pairs of distinct locations of the warehouse. The dynamics of the problem is modelled through a space-time network ${\mathcal {G}}=({\mathcal {N}}, {\mathcal {A}})$. Specifically, the time horizon is discretized into T time periods of equal length through $T+1$ time instants. The set ${\mathcal {N}}^P$ is then replicated $T+1$ times, resulting in set ${\mathcal {N}}$. A node in ${\mathcal {N}}$ is defined by a couple (i, t), with $i \in {\mathcal {N}}^P$ and $t \in \{0,\dots ,T\}$, and represents one of the locations of the warehouse at one of the considered $T+1$ time instants. The set ${\mathcal {A}}$ is composed of two subsets: the subset ${\mathcal {A}}^{H}$ of holding arcs, i.e., arcs of form $((i,t), (i,t+1))$, for any $i \in {\mathcal {N}}^P$ and $t \in \{0, \dots , T-1\}$, used to model idle time of items or vehicles in a given node for one time period, and the subset ${\mathcal {A}}^{M}$ of moving arcs, composed of arcs of form $((i,t), (j,t'))$ with $(i,j) \in {\mathcal {A}}^P$, $t \in \{0, \dots , T-\tau _{i,j}\}$ and $t'=t+\tau _{i,j}$, where $\tau _{i,j}$ denotes the travel time from i to j in the physical network ${\mathcal {G}}^P$. ${\mathcal {A}}^{M}$ thus models movements of items or vehicles between two different locations in different time periods.

4.2 Battery management

We defined an approximate energy consumption model to compute the energy consumption of an EV within the warehouse. This energy consumption model has been used in several scientific papers to estimate the energy consumption of electric trucks, electric cars, and electric buses (Basso et al. 2019b; Macrina et al. 2019; Travesset-Baro et al. 2015; Tang et al. 2019; Gao et al. 2017). Multiple factors may affect the energy consumption rate. However, loading items and moving across the warehouse can be regarded as the two most important sources of energy in warehousing activities.

The energy consumption for movements within the warehouse is computed in three steps, as in Fiori et al. (2016). In step one, the mechanical power $P_W$ needed at the wheels to make the vehicle move or brake, also called power at the wheels, is computed. In step two, the electric power $P_E$ needed to achieve the mechanical power $P_W$ is calculated. Finally, in step three, the required electric power $P_E$ is converted into the power $P_B$ that has to be taken from the battery.

To calculate the mechanical power, we use the well-known equation (1), based on Newton’s Law of Motion applied to move vehicles, which has been utilized by several authors in different logistic contexts. It computes, as a function of time, the mechanical power needed to overcome the force of air resistance, the rolling resistance, and the force of gravity acting against the motion of the vehicle:

$$\begin{aligned} P_{W}(t) = \left[ m\,a(t) + m\,g\, cos(\theta ) \, C_r + \frac{1}{2}\rho _{Air} \, A_f \, C_D \, v^2(t) + m \, g \, sin(\theta ) \right] \, v(t), \end{aligned}$$

(1)

where $m=m_0+w$ is the vehicle mass [kg], which includes the curb weight $m_0$ of the vehicle and the weight w of the transported items, a(t) is the acceleration of the vehicle [$m/s^2$], g is the gravitational acceleration constant [$m/s^2$], $\theta$ is the road inclination angle, $C_r$ is the rolling resistance parameter [DN] which is a function of the road surface type, road condition, and vehicle tire type, $\rho _{Air}$ is the density of air [$kg/m^3$], $A_f$ is the frontal area of the vehicle [$m^2$], $C_D$ is the aerodynamic drag resistance coefficient of the vehicle [DN], and v(t) is the vehicle speed [m/s]. When (1) is computed second-by-second, $P_{W}(t)$ may be expressed in Joules.

Given the mechanical power, the amount of electric power that the electric motor needs to produce has then to be found. The relationship $P_{E}(t)=P_W(t)/\eta$ can be used as an approximation, where $\eta$ represents the efficiency of the electric motor in limiting the losses in energy conversion. Nevertheless, the vehicles we consider are equipped with the regenerative braking recovery mechanism, that allows converting the kinetic energy produced when the vehicle slows down into electrical energy to recharge the battery. The negative sign of $P_W(t)$ (due to the sign of a(t) or due to the road grade) determines when energy is recovered, whose recovered amount depends on the efficiency of the recovery mechanism of the electric motor. Therefore, the above-mentioned relationship can be extended as follows:

$$\begin{aligned} P_{E}(t)= {\left\{ \begin{array}{ll} P_W(t)/\eta ^{+} & \quad {}\hbox {if } P_{W}(t) \ge 0, \\ P_W(t)/\eta ^{-} & \quad {}\hbox {if }P_{W}(t) < 0. \end{array}\right. } \end{aligned}$$

(2)

This gives both the amount of electric power needed by the electric motor to provide the mechanical power required to move the vehicle, given its efficiency in converting electrical power into mechanical power (case $\eta ^{+} \le 1$), and also the amount of recovered energy when the vehicle brakes, given the efficiency of the recovery mechanism (case $\eta ^{-} \ge 1$). Thanks to the second-by-second calculation, $P_{E}(t)$ may be still expressed in Joules.

Finally, we approximate the battery power as follows:

$$\begin{aligned} P_{B}(t)= {\left\{ \begin{array}{ll} P_E(t)/\phi ^{+} & \quad {}\hbox {if } P_{E}(t) \ge 0, \\ P_E(t)/\phi ^{-} & \quad {}\hbox {if }P_{E}(t) < 0. \end{array}\right. } \end{aligned}$$

(3)

When $P_{E} \ge 0$, then $P_{B}$ represents the battery power that is necessary to obtain the required electric energy $P_{E}$ to move the vehicle, thus the battery is discharged. When $P_{E} < 0$ instead, $P_{B}$ models the recovery of energy that is produced when the vehicle brakes, thus the battery is charged. The battery has losses during both charging and discharging and $\phi ^{+}$ and $\phi ^{-}$ thus represent the battery efficiency. Also $P_{B}(t)$ can be expressed in Joules.

To calculate the battery power $P^L_{B}$ needed to lift items from the ground, we first consider the mechanical power needed to lift the object through the well-known equation

$$\begin{aligned} P^L_{W} = w \, g \, h, \end{aligned}$$

where w is the mass [kg] of the item to lift, g is the gravitational acceleration constant [$m/s^2$] and h is the height to which the item is lifted [m]. We then approximate $P^L_{B}$, required to perform this operation, by following the same reasoning above described and considering the energy losses that occur when converting mechanical power into electric power and then into battery power:

$$\begin{aligned} P^L_{B}= \frac{P^L_W}{\eta ^{+}\phi ^{+}}. \end{aligned}$$

(4)

On the other hand, the mechanical power required to lower items can be considered negligible.

To take into account the battery energy consumption in the model formulation, we define two parameters for each arc $(i,j)\in {\mathcal {A}}^{P}$, which are related to the discharge of the battery of an EV when moving along the arc. The parameter $e^{ij}$ represents the battery energy consumed by the vehicle to move empty along (i, j), while the parameter $e^{ijk}$, with $k \in {\mathcal {K}}$, represents the additional battery energy consumed by the vehicle to move along (i, j) if it is loaded with items of product type k, per unit of load. More in detail, to assign a value to $e^{ij}$, we first split the associated travel time $\tau _{ij}$ into three segments, as in Fig. 2. In the first segment, from 0 to $\tau '$, the vehicle is supposed to accelerate till its maximum allowed speed, $v_{max}$, thus a positive value for a(t) is considered in (1); in the second segment, from $\tau '$ to $\tau ''$, the vehicle is assumed to travel at a constant speed, thus $a(t) = 0$ in (1); in the third segment, from $\tau ''$ to $\tau _{ij}$, the vehicle brakes and a(t) has a negative value in (1). Given the application context under consideration, we assume that a vehicle can always reach its maximum speed when travelling along an arc (see Macrina et al. 2019, for an application in which vehicles may not always reach their maximum speed on a segment). The battery power associated with (i, j) is defined by the sum of the battery powers required by the vehicle for each of the three segments above described. The amount of battery power associated with each segment is computed through formulas (1), (2) and (3), by integrating over the corresponding length of time (i.e, for the first segment the integration is from 0 to $\tau '$, for the second segment from $\tau '$ to $\tau ''$ and the third segment from $\tau ''$ to $\tau _{ij}$). Notice that $e^{ij}$ refers to the battery energy consumed by an empty vehicle, thus $m = m_0$ in formula (1). The same kind of procedure has been followed to compute the value of parameter $e^{ijk}$ for $(i,j)\in {\mathcal {A}}^{P}$ and $k\in {\mathcal {K}}$, now setting $m = w$ in (1), with w denoting the weight of an item of type k, in order to determine the additional energy required to move with a load of type k.

Fig. 2

Speed profile example

Moreover, we define further parameters, i.e., $e^k$ for each $k \in {\mathcal {K}}$, to model the energy consumed by an EV to lift one unit of product type k. Parameters $e^k$ are computed through formula (4). Notice that lifting operations can be performed only at nodes in $S^k_{in} \cup {\mathcal {S}}^k_{out} \cup {\mathcal {B}}$.

Regarding the operations at the charging station, a parameter $e^r$ is associated with the holding arcs connecting node c, which represents the charging station, at two consecutive time instants, denoting the increase of the battery energy, for one period of time, if the vehicle recharges at c. A maximum number $c_{max}$ of EVs can be recharged at the same time. In addition, we denote with $[B^{-}, B^{+}]$ the range in which the charge of each battery must always be maintained, with $z^v_0\in [B^{-}, B^{+}]$ representing the charge of vehicle $v \in {\mathcal {V}}^{E}$ at the beginning of the time horizon.

Finally, two additional parameters, $\Delta \in (0,1)$ and $\Theta \in [B^{-}, B^{+}]$, are introduced to regulate the state of charge of the EVs at the end of the time horizon. Specifically, they are introduced to guarantee that the fleet of EVs must not be completely discharged at the end of the considered time horizon, but rather must have an amount of charge adequate to be operative in the next time horizon. Application contexts where such a requirement is extremely relevant and useful are warehouses performing all-day-long operations (i.e., 24 h per day), like in the case study addressed. As detailed in Sect. 4.3, this behaviour is guaranteed by imposing that the cumulative state of charge of the batteries of all the EVs at the end of the time horizon be not lower than a percentage of the initial total charge of the fleet. Such a percentage is defined by the parameter $\Delta$. This is quite flexible, indeed, since the level of charge of any EV at the end of the time horizon can be suitably chosen by the solver, provided that the minimum total charge for the fleet is guaranteed. Nevertheless, a minimum charge for each EV at the end of the time horizon is imposed too, through parameter $\Theta$. Specifically, such a parameter is introduced to ensure that each EV, at the beginning of the next time horizon, has enough charge to move to the charging station, thus preventing its block.

4.3 Energy constraints

GSRP is formulated in terms of multicommodity flows on the space-time network ${\mathcal {G}}$, and a Mixed-Integer Linear Programming (MILP) model on ${\mathcal {G}}$ is proposed.

Let us introduce the main families of variables used. The following four families of variables model the routing of vehicles and commodities along the network. In the variable definition, ${\mathcal {A}}_{F1}$, ${\mathcal {A}}_{F2}$, ${\mathcal {A}}_{in}$ and ${\mathcal {A}}_{out}$ denote the subsets of arcs where vehicles of type F1, vehicles of type F2, incoming items and outgoing items can move, respectively:

• $x^v_{(i,t)(j,t')} \in \lbrace 0,1 \rbrace$, for any $v \in {\mathcal {V}}^{1}$ and $((i,t),(j,t')) \in {\mathcal {A}}_{F1}$, indicates whether the vehicle v of type F1 passes on the arc $((i,t),(j,t'))$ or not;

• $x^v_{(i,t)(j,t')} \in \lbrace 0,1 \rbrace$, for any $v \in {\mathcal {V}}^{2}$ and $((i,t),(j,t')) \in {\mathcal {A}}_{F2}$, indicates whether the vehicle v of type F2 passes on the arc $((i,t),(j,t'))$ or not;

• $y^k_{(i,t)(j,t')} \in \mathbb {Z}_+$, for any $k \in {\mathcal {K}}_{in}$ and $((i,t),(j,t')) \in {\mathcal {A}}_{in}$, indicates the number of incoming items of product type k passing on the arc $((i,t),(j,t'))$;

• $y^k_{(i,t)(j,t')} \in \mathbb {Z}_+$, for any $k \in {\mathcal {K}}_{out}$ and $((i,t),(j,t')) \in {\mathcal {A}}_{out}$, indicates the number of outgoing items of product type k passing on the arc $((i,t),(j,t'))$.

Moreover, regarding battery management, we define:

• $z^v_t \in \mathbb {R}_+$, for any $v \in {\mathcal {V}}^{E}$ and $t \in \{1, \dots , T\}$, indicates the state of charge of the battery of the EV v at time t.

The objective function of the MILP model is defined as follows:

$$\begin{aligned} \begin{array}{rl} \min \ {} &{} \displaystyle \sum _{v \in {\mathcal {V}}^{1}} \sum _{\begin{array}{c} ((i,t),(j,t')) \in {\mathcal {A}}_{F1}:\\ i\ne \omega ^{1},\,j \ne \omega ^{1} \end{array}} \tau _{i,j} \, x^v_{(i,t)(j,t')} + \sum _{v \in {\mathcal {V}}^{2}} \sum _{\begin{array}{c} ((i,t),(j,t')) \in {\mathcal {A}}_{F2}:\\ i\ne \omega ^{2},\,j \ne \omega ^{2} \\ i,\,j \ne c \end{array}} \tau _{i,j} \, x^v_{(i,t)(j,t')}\\ &{} \displaystyle + \psi \sum _{k \in {\mathcal {K}}_{in}} \sum _{\begin{array}{c} ((i,t),(j,t')) \in {\mathcal {A}}_{in}: \\ i,\,j \in {\mathcal {R}} \end{array}} y^k_{(i,t)(j,t')} + \xi \sum _{k \in {\mathcal {K}}_{out}} P^k + \zeta \sum _{v \in {\mathcal {V}}^{E}} \sum _{t=0}^{T-1} x^v_{(c,t)(c,t+1)}. \end{array} \end{aligned}$$

(5)

It is composed of five summations. The first two summations define the primary optimization goal, i.e., minimizing the travel time of all the vehicles in ${\mathcal {V}}^{1}$ and of all the vehicles in ${\mathcal {V}}^{2}$ within the warehouse. Arcs entering or leaving the parking areas are not considered to encourage vehicles to come back to their parking areas when idle, so limiting congestion along the network. The other three terms are used to handle additional features to consider during the optimization process, which are thus modelled through soft objectives. Specifically, the third summation aims at penalizing the time of permanence of the incoming items on the input points and, consequently, to favour the movements of such items towards other spots of the warehouse. The anticipation movement of outgoing items toward the collection area is addressed through the fourth summation. Being ${\mathcal {N}}^-(\pi , t')$ the set of nodes linked to $(\pi ,t')$ via an entering arc, the penalty $P^k$, for any $k \in {\mathcal {K}}_{out}$, is defined as follows:

$$\begin{aligned} P^k = \max \left\{ 0 , \ \sum _{t=0}^{\tilde{T}}d_{out}^k(\pi ,t) - \left[ u^k_{\pi } + \sum _{t=0}^{T}\sum _{(j,t) \in {\mathcal {N}}^{-}(\pi ,t')} y^k_{(j,t)(\pi ,t')} \right] \right\} , \end{aligned}$$

(6)

where $u^k_{\pi }$ is the number of items of product type k in the collection area at the beginning of the time horizon, $d_{out}^k(\pi ,t)$ is the number of items of type k which are requested in the collection area at the latest time t, and ${\tilde{T}} > T$ is an extended time instant to be considered for the anticipation moves. Therefore, if during the considered time horizon the number of outgoing items of type k moved to the collection area is not enough to satisfy both the demand of k in the time horizon and also in the extended one, then the penalty is set proportionally to the unmet demand; otherwise, the penalty is equal to 0. Finally, the fifth summation addresses the issue of unnecessary recharging operations, i.e., when an EV is charged even though it does not have any forthcoming operation planned to perform. Unnecessary recharge periods may raise the cost of warehouse management, uselessly consume resources within the warehouse, and increase the environmental impact (Bartolini et al. 2019). The three soft objectives are weighted through parameters $\psi$, $\xi$, and $\zeta$, respectively, to state their mutual priorities and their priority with respect to the primary optimization goal, and also to allow a comparison between the different units of measure of the soft and primary criteria, through a suitable parameter calibration (for more details see Lanza et al. 2022d).

Since peculiar to GSRP, we report below only those constraints which regulate the energy behaviour of the battery of an EV plus some additional energy constraints. The constraints in the first group are presented separately for the case of level-of-charge decrease (when the vehicle travels along a moving arc as well as when the vehicle lifts some items), level-of-charge increase (when the vehicle idles on a holding arc corresponding to the charging station), and constant level-of-charge (when the vehicle idles in a location other than the charging station). In such constraints, M is an input parameter defined as $M=\,B^{+}-B^{-}$. Moreover, $L_j$ is a parameter equal to 1 if $j \in S^k_{in} \cup {\mathcal {S}}^k_{out} \cup {\mathcal {B}}$, and 0 otherwise.

Battery level-of-charge constraints

$$\begin{aligned} \displaystyle z^v_{t'} \le\,&\displaystyle z^v_{t} - e^{ij} x^v_{(i,t)(j,t')} - \sum _{k \in \mathcal {K}} e^{ijk} y^k_{(i,t)(j,t')} - L_j \sum _{k \in \mathcal {K}} e^{k} y^k_{(i,t)(j,t')}\nonumber \\&\quad + M\left[ 1-x^v_{(i,t)(j,t')}\right] \qquad \forall \ v \in \mathcal {V}^E,\ \forall \ ((i,t)(j,t')) \in \mathcal {A}_{F2}:\ i\ne j, \end{aligned}$$

(7)

$$\begin{aligned} z^v_{t'} \ge\,&\displaystyle z^v_{t} - e^{ij} x^v_{(i,t)(j,t')} - \sum _{k \in \mathcal {K}} e^{ijk} y^k_{(i,t)(j,t')} - L_j \sum _{k \in \mathcal {K}} e^{k} y^k_{(i,t)(j,t')}\nonumber \\&\quad - M\left[ 1-x^v_{(i,t)(j,t')}\right] \qquad \forall \ v \in \mathcal {V}^E,\ \forall \ ((i,t)(j,t')) \in \mathcal {A}_{F2}:\ i\ne j. \end{aligned}$$

(8)

Constraints (7)–(8) model the discharge of the battery. Recalling that a moving arc can be used by at most one vehicle at a time, when a vehicle v travels along a moving arc $((i,t)(j,t'))$ these constraints imply

$$\begin{aligned} z^v_{t'} = z^v_{t} - e^{ij} - \sum _{k \in \mathcal {K}} e^{ijk} y^k_{(i,t)(j,t')} - L_j \sum _{k \in \mathcal {K}} e^{k} y^k_{(i,t)(j,t')}, \end{aligned}$$

thus defining the level of charge of its battery at the time instant $t'$ as the level of charge of its battery at the time instant t minus the energy necessary to v to move empty on the arc (the second term in the equation), the additional energy used if the vehicle is loaded (the third term in the equation) and the energy used to lift items at location j, if necessary (the last term in the equation). If the vehicle does not travel along the arc, then constraints (7)–(8) are satisfied since weaker than constraints (15).

$$\begin{aligned}&z^v_{t+1} \le z^v_{t} + e^r x^v_{(c,t)(c,t+1)} + M\left[ 1-x^v_{(c,t)(c,t+1)}\right] \qquad \forall \ v \in \mathcal {V}^E,\ \forall \ ((c,t)(c,t+1)) \in \mathcal {A}_{F2}, \end{aligned}$$

(9)

$$\begin{aligned}&z^v_{t+1} \ge z^v_{t} + e^r x^v_{(c,t)(c,t+1)} - M\left[ 1-x^v_{(c,t)(c,t+1)}\right] \qquad \forall \ v \in \mathcal {V}^E,\ \forall \ ((c,t)(c,t+1)) \in \mathcal {A}_{F2}. \end{aligned}$$

(10)

Constraints (9)–(10) model the recharge of the battery. When the vehicle is at the charging station c, i.e., it is on a holding arc of form $((c,t)(c,t+1))$, these constraints define the level of charge of the battery at time instant $t+1$ as the level of charge of the battery at the time instant t plus the energy recharged during one time period at c.

$$\begin{aligned}&z^v_{t+1} \le z^v_{t} + M\left[ 1-x^v_{(i,t)(i,t+1)}\right] \qquad \forall \ v \in \mathcal {V}^E,\ \forall \ ((i,t)(i,t+1)) \in \mathcal {A}_{F2}:\ i\ne c, \end{aligned}$$

(11)

$$\begin{aligned}&z^v_{t+1} \ge z^v_{t} - M\left[ 1-x^v_{(i,t)(i,t+1)}\right] \qquad \forall \ v \in \mathcal {V}^E,\ \forall \ ((i,t)(i,t+1)) \in \mathcal {A}_{F2}:\ i\ne c. \end{aligned}$$

(12)

Finally, constraints (11)–(12) state that the level of charge of the battery is unchanged if the vehicle is idling on a location of the warehouse other than c.

Additional energy constraints

$$\begin{aligned}&\sum _{v \in \mathcal {V}^E} x^v_{(c,t)(c,t+1)} \le c_{max}&\qquad \forall \ t \ge 0, \end{aligned}$$

(13)

$$\begin{aligned} \sum _{v \in \mathcal {V}^E} z^v_{T} \ge \Delta \sum _{v \in \mathcal {V}^E}z^v_0, \end{aligned}$$

(14)

$$\begin{aligned}&B^{-} \le z^v_t \le B^{+}&\qquad \forall \ v \in \mathcal {V}^E, \ \forall \ t \ge 1, \end{aligned}$$

(15)

$$\begin{aligned}&z^v_T \ge \Theta&\qquad \forall \ v \in {\mathcal {V}}^{E}. \end{aligned}$$

(16)

Constraints (13) define the maximum number of vehicles that may be charged simultaneously at the charging station. Constraint (14) imposes that, at the end of the time horizon, the total charge of all the EVs must be greater than or equal to a prespecified percentage $\Delta$ of the initial total battery charge, i.e., $\Delta \sum _{v \in \mathcal {V}^E}z^v_0$. This requirement is particularly useful when the fleet of EVs must be prompt to perform operations in the next time horizon, for instance in applications where operations must be performed 24 h per day. Constraints (15) define the lower and upper thresholds for the ideal operating conditions of the battery. Finally, constraints (16) impose that the level of charge of each EV at the end of the time horizon be greater than or equal to a minimum threshold $\Theta$. Notice that, when a unique EV is considered and $\Delta > \Theta$, then constraints (16) are indeed redundant. However, when several EVs are used, constraints (16) ensure that, at the beginning of the next time horizon, each vehicle has a state of charge sufficient to perform some basic operations rather than being completely discharged. As observed in Sect. 1, the additional energy constraints and the last term of the objective function (5) render the studied problem still more adherent to the reality of warehouses such as the one analysed in the considered case study.

The mathematical formulation is completed by flow conservation constraints for the incoming and the outgoing items as well as for the vehicles, by constraints ensuring their correct transportation and routing within the warehouse, by linking capacity constraints for vehicles and incoming and outgoing flows, by constraints ensuring the respect of the due dates for the outgoing items, by location capacity constraints, and by constraints related to some of the soft objectives pursued, i.e., the application of the strict management policy for picking and put-away operations, and the security requirements for workers. Such additional constraints, presented in Lanza et al. (2022d) and referring to a traditional context with only conventional vehicles, are reported in Appendix A for the sake of completeness. We still emphasize that the proposed formulation extends the one in Lanza et al. (2022c) since additional battery management features (constraints (14)), the capacity constraints at the charging station (constraints (13)), and a further soft objective (given by the fifth term in (5)) are considered.

5 Resolution approaches

Even if the problem is NP-hard, the state-of-the-art commercial solver CPLEX was able to determine an optimal solution, or a suitable lower bound, to small to medium size instances, as reported in Sect. 6.4.1. On the other hand, the commercial solver applied to the mathematical model was not able to solve larger instances, due to the problem complexity. In particular, for medium to large-scale instances, like the real ones, the proposed formulation may have a very high dimension because of the huge number of products and storage locations involved in storing and retrieving operations, and the constraints related to the battery management. This has motivated the proposal of matheuristic approaches to efficiently determine solutions of good quality also in case of instances of larger size.

In this section, we describe the matheuristic approaches we designed to solve GSRP, whose performance is evaluated in the experimental sections. More precisely, in Sect. 5.1 we present the main features of the time decomposition matheuristic proposed in Lanza et al. (2022c) and then suitably extended to the problem addressed, in Sect. 5.2 we describe a new approach based on activity decomposition, while in Sect. 5.3 we propose a further approach, obtained by a suitable combination of the two before mentioned approaches, which proved to be particularly efficient on the tested instances.

5.1 Time horizon decomposition matheuristic

The time horizon decomposition approach has been initially proposed for SRP, which showed a very good performance, as detailed in Lanza et al. (2022d). Then, it has been developed a modified version of the algorithm to deal with the GSRP in Lanza et al. (2022c). Here, the time horizon decomposition approach described in Lanza et al. (2022c) has been further extended to take into account the new green features of the problem currently addressed.

In this approach, the planning horizon is divided into $\Lambda$ subperiods of equal length. Each subperiod gives rise to a subproblem, whose features are those of the original problem restricted to the considered subperiod. The $\Lambda$ subproblems are then sequentially solved by using a MILP solver in such a way that the final state of the system obtained solving subproblem $\lambda -1$ becomes the initial state of the system when solving subproblem $\lambda$, for $\lambda = 2,\dots ,\Lambda$. Specifically, the state of the system in each subproblem takes into account the position of vehicles and items within the warehouse and the level of charge of each battery. In particular, the level of charge of an EV v at the initial time instant $0_{\lambda }$ of subproblem $\lambda$ is defined as $z^v_{0_{\lambda }} = z^v_{T_{\lambda -1}}$, where $z^v_{T_{\lambda -1}}$ is the level of charge of the battery of v at the final time instant $T_{\lambda -1}$ of subproblem $\lambda -1$.

Once the $\Lambda$ subproblems are solved, to construct a solution for GSRP it is sufficient to concatenate the $\Lambda$ solutions in an increasing order with respect to the subperiod addressed, i.e., from subperiod 1 to subperiod $\Lambda$. The matheuristic approach, hereafter referred to as TD, is summarized in Algorithm TD.

5.2 Activity decomposition matheuristic

The activity decomposition matheuristic consists of two phases. Firstly, a subproblem restricted to the incoming items, called $P_{K_{in}}$, is solved to define the complete itineraries (i.e., from the conveyors towards the preassigned storage locations over the complete time horizon) for each $k \in {\mathcal {K}}_{in}$, as well as to define the complete routing for each vehicle $v \in {\mathcal {V}}^{1}$. This subproblem also determines a partial routing for each vehicle $v \in {\mathcal {V}}^{2}$, only concerning the operations which are required to move the products in ${\mathcal {K}}_{in}$ within the storage area. Then, a subproblem restricted to the outgoing items, called $P_{K_{out}}$, is solved to define the complete itineraries (i.e., from the storage locations towards the collection area over the complete time horizon) for each $k \in {\mathcal {K}}_{out}$, as well as the complementary routing for each $v \in {\mathcal {V}}^{2}$, only concerning the operations required to move the products in ${\mathcal {K}}_{out}$. The two subproblems are formulated on two reduced auxiliary graphs, built starting from ${\mathcal {G}}=({\mathcal {N}}, {\mathcal {A}})$, and solved sequentially by fixing, in $P_{K_{out}}$, some of the variables related to the routing of the vehicles in ${\mathcal {V}}^{2}$ when performing put-away operations, as given by the resolution of $P_{K_{in}}$. The solution to GSRP is then given by the routing of the vehicles in ${\mathcal {V}}^{1}$ and the itinerary of the products in $K_{in}$ obtained by the first phase, by the itinerary of the products in $K_{out}$ returned by the second phase, and by the routing of the vehicles in ${\mathcal {V}}^{2}$ obtained by the combinations of the two partial routings returned by $P_{K_{in}}$ and $P_{K_{out}}$, respectively.

The activity decomposition matheuristic, referred to as AD, is summarized in Algorithm AD. The two phases are detailed in the following.

5.2.1 Subproblem $P_{K_{in}}$

$P_{K_{in}}$ is defined over a reduced auxiliary graph, $\tilde{{\mathcal {G}}}_{K_{in}}=(\tilde{{\mathcal {N}}}_{K_{in}}, \tilde{{\mathcal {A}}}_{K_{in}})$, which is obtained from the complete auxiliary graph ${\mathcal {G}}=({\mathcal {N}}, {\mathcal {A}})$ by removing all the nodes in the set ${\mathcal {S}}^k_{out}$ and node $\pi$, as well as all the arcs incident to such nodes. Therefore, $\tilde{{\mathcal {N}}}_{K_{in}} = {\mathcal {R}} \cup {\mathcal {B}} \cup \mathcal {S}^k_{in} \cup \{\omega ^{1} \} \cup \{ \omega ^{2} \} \cup \{ c \}$, while $\tilde{{\mathcal {A}}}_{K_{in}}$ is composed of all the arcs $((i,t)(j,t')) \in {\mathcal {A}}$ whose end nodes belong to $\tilde{{\mathcal {N}}}_{K_{in}}$. In other words, $\tilde{{\mathcal {A}}}_{K_{in}}$ contains only the arcs which are necessary for the transportation of the product types in ${\mathcal {K}}_{in}$ from the conveyors towards their preassigned storage locations.

$P_{K_{in}}$ is mathematically formulated by using a subset of the variables introduced in Sect. 4. Specifically, we only consider variables $x^v_{(i,t)(j,t')} \in \lbrace 0,1 \rbrace$, for any $v \in {\mathcal {V}}^{1}$ and $((i,t),(j,t')) \in {\mathcal {A}}_{F1}$, variables $x^v_{(i,t)(j,t')} \in \lbrace 0,1 \rbrace$, for any $v \in {\mathcal {V}}^{2}$ and $((i,t),(j,t')) \in {\mathcal {A}}_{F2} \cap \tilde{{\mathcal {A}}}_{K_{in}}$, and variables $y^k_{(i,t)(j,t')} \in \mathbb {Z}_+$, for any $k \in {\mathcal {K}}_{in}$ and $((i,t),(j,t')) \in {\mathcal {A}}_{in}$, plus some auxiliary variables related to the management policy for the put-away operations (i.e., the variables $\alpha (s^k, t) \in \lbrace 0,1 \rbrace$, for any $s^k \in {\mathcal {S}}^k_{in}$, $k \in {\mathcal {K}}_{in}$, and $t=0,\dots ,T$ introduced in Appendix A). The $P_{K_{in}}$ formulation thus consists of the objective function (5), without terms (6), and by constraints (7)–(16), (18)–(21), (22), (26)–(29), (31), (33)–(35), where for constraints from (18) to (35) we refer to the Appendix A. Notice that such constraints must indeed be restricted to the subset of variables mentioned before.

The solution of $P_{K_{in}}$ thus provides the complete schedule of the operations of the fleet F1 and the partial schedule of the operations of the fleet F2 involving only storage operations, plus the complete itinerary of the incoming items towards their storage locations.

5.2.2 Subproblem $P_{K_{out}}$

Also $P_{K_{out}}$ is defined over a reduced auxiliary graph, $\tilde{{\mathcal {G}}}_{K_{out}}=(\tilde{{\mathcal {N}}}_{K_{out}}, \tilde{{\mathcal {A}}}_{K_{out}})$, which is obtained from the complete auxiliary graph ${\mathcal {G}}=({\mathcal {N}}, {\mathcal {A}})$ by removing the nodes in the set ${\mathcal {R}}$, the parking node $\omega ^{1}$, as well as all the arcs incident to such nodes. Specifically, $\tilde{{\mathcal {N}}}_{K_{out}} = {\mathcal {B}} \cup \mathcal {S}^k_{in} \cup {\mathcal {S}}^k_{out} \cup \{\pi \} \cup \{ \omega ^{2} \} \cup \{c\}$, while the set $\tilde{{\mathcal {A}}}_{K_{out}}$ is composed of all the arcs $((i,t)(j,t')) \in {\mathcal {A}}$ whose end nodes belong to $\tilde{{\mathcal {N}}}_{K_{out}}$. The set thus includes all the arcs which are necessary to move the product types in ${\mathcal {K}}_{out}$ from their storage locations towards the collection area, and the subset of arcs used by the vehicles of type F2 to perform picking operations.

Also $P_{K_{out}}$ is mathematically formulated by means of a subset of the variables introduced in Sect. 4. Specifically, we only consider variables $x^v_{(i,t)(j,t')} \in \lbrace 0,1 \rbrace$, for any $v \in {\mathcal {V}}^{2}$ and $((i,t),(j,t')) \in {\mathcal {A}}_{F2}$, variables $y^k_{(i,t)(j,t')} \in \mathbb {Z}_+$, for any $k \in {\mathcal {K}}_{out}$ and $((i,t),(j,t')) \in {\mathcal {A}}_{out}$, plus some auxiliary variables related to the management policy for the picking operations (i.e., the variables $\beta (s^k, t) \in \lbrace 0,1 \rbrace$, for any $s^k \in {\mathcal {S}}^k_{out}$, $k \in {\mathcal {K}}_{out}$, and $t = 0, \ldots , T$ introduced in Appendix A.) The $P_{K_{out}}$ formulation thus consists of the objective function (5), without the term related to the idle time of incoming items on the conveyors (i.e., the third one), and constraints (7)–(16), (19), (21), (23)–(25), (27), (29), (30), (32) and (36)–(38), where for constraints from (19) to (38) we refer again to Appendix A. Indeed, in order to incorporate the partial schedule of the vehicles in F2 returned by the first phase of the approach, in the $P_{K_{out}}$ formulation we fix to 1 the value of the variables $x^v_{(i,t)(j,t')}$ corresponding to the arcs travelled by a vehicle of type F2 when engaged in some storing operations, and we fix all the integer variables $y^k_{(i,t)(j,t')}$, for $k \in {\mathcal {K}}_{in}$, as returned by the solution of subproblem $P_{K_{in}}$. The solution of $P_{K_{out}}$ thus provides the complete schedule of the operations of the fleet F2, involving both put-away and picking operations, as well as the complete itinerary for outgoing items towards the collection area.

5.3 Activity and time horizon decomposition matheuristic

This approach combines the features of the two algorithms presented before, in a nested way. The original planning horizon is divided into $\Lambda$ subperiods of equal length, by applying the time horizon decomposition. Each subperiod, say $\lambda$, gives rise to a temporal subproblem, whose features are those of the original problem restricted to the considered subperiod. The resolution of each temporal subproblem is then performed in two phases, by applying the activity decomposition. Specifically, in the first phase of the temporal subproblem resolution, an activity subproblem $P^{\lambda }_{K_{in}}$, related to the temporal subproblem taken into account, is solved. Then, an activity subproblem $P^{\lambda }_{K_{out}}$, always related to the temporal subproblem under consideration, is solved by fixing to 1 the value of the variables $x^v_{(i,t)(j,t')}$ corresponding to the arcs travelled by a vehicle of type F2 when performing storing operations, and fixing all the integer variables $y^k_{(i,t)(j,t')}$, for $k \in {\mathcal {K}}_{in}$, as returned by the resolution of $P^{\lambda }_{K_{in}}$.

The constraints modelling $P^{\lambda }_{K_{in}}$ and $P^{\lambda }_{K_{out}}$ are the same as in AD, restricted however to the considered subperiod $\lambda$, with the only exception of constraint (14), which is included also in the formulation related to the subperiods $\lambda = 1,\dots , \Lambda -1$, although in a modified form. Recall that (14) imposes that, at the end of the overall time horizon, the total charge of the batteries must be not lower than a prespecified level. Including such a constraint only in the last subproblem, i.e., in the one related to subperiod $\Lambda$, may generate very hard, or even infeasible, subproblems, for example when the time required to charge vehicles till the required level $\Delta \sum _{v \in \mathcal {V}^E}z^v_0$ is longer than the length of subproblem $\Lambda$. For this reason, the following version of constraint (14) has been included in the subproblems related to $\lambda = 1,\dots ,\Lambda -1$:

$$\begin{aligned} \sum _{v \in {V}^E} z^v_{T_{\lambda }} \ge \frac{\lambda \Delta }{\Lambda } \sum _{v \in \mathcal {V}^E}z^v_0, \end{aligned}$$

(17)

where $T_{\lambda }$ denotes the final time instant of subperiod $\lambda$. Constraint (17) imposes that a fraction of the required final level of charge is satisfied at the end of each subperiod. This allows one not to concentrate the charging operations in the last subperiod, but rather to split them among the subperiods.

The $\Lambda$ temporal subproblems are sequentially solved in such a way that the final state of the system when solving subproblem $\lambda -1$ becomes the initial state of the system when solving subproblem $\lambda$, for $\lambda = 2,\dots ,\Lambda$. In order to construct a solution for GSRP, it is then necessary to concatenate the $\Lambda$ solutions in increasing order with respect to the subperiod addressed, i.e., from subperiod 1 to subperiod $\Lambda$, as in TD. Hereafter we refer to this matheuristic as ATD, and we summarize it in Algorithm ATD.

We conclude by emphasizing that algorithms AD and ATD are indeed different, although both are based on the common idea to suitably decompose the mathematical model, in order to solve efficiently also instances of large size. The algorithm AD, in fact, rather than jointly addressing put-away and picking operations, like in the mathematical model, firstly takes decisions about the put-away activities - the most critical ones, given the policy adopted - and after about the picking activities. This is done by looking at the overall planning period. The algorithm ATD instead does not consider the entire planning horizon to schedule the two types of operations, rather it decomposes the planning horizon into subperiods, and then applies the philosophy of AD to each subperiod, by managing the different subperiods in a chronological order. That is, AD performs only an activity decomposition, while ATD performs two nested levels of decompositions: a time decomposition and then, for each subperiod, an activity decomposition. The two matheuristics have thus different characteristics. Thanks to the double decomposition, the dimension of the subproblems of ATD is smaller than in AD, and this may facilitate the problem resolution. However, ATD has a more restricted vision of the future, due to the performed time decomposition, and therefore in principle it could return solutions of worse quality compared to AD. Instead, the computational analysis reported in Sects. 6.4 and 6.5 has shown that ATD is able to efficiently combine the temporal and the activity decomposition by computing solutions of good quality to GSRP.

6 Numerical experiments

6.1 The case study

The case study refers to a production site related to the tissue sector located in Tuscany. The site comprises a production area, a warehouse larger than 10,000 $\hbox {m}^2$, a collection area, and several shipping docks. In turn, the warehouse is composed of 4 departments. Each department has a rectangular internal layout with a certain number of parallel narrow storage aisles and parallel wide cross aisles. The storage area is thus divided into blocks of storage locations framed by the aisles. Items are homogeneously stored back-to-back in each storage location, thus generating horizontal stacks accessible only frontally. There are 858 stacks organized into 29 blocks. Following the pick-and-sort policy, the collection area is used to gather picked items and establish order integrity before loading trucks. It can stock up to 700 items and is normally filled up as much as possible during the night to quickly start the truck-loading operations the next morning. The warehouse also contains 6 collectors at the entrance of the departments, whose capacities range from 2 to 8 items.

The production site works daily on three shifts of 8 h, by producing more than 300 different types of tissue products for sanitary and domestic use. Put-away operations are performed 24 h non-stop in the warehouse. Items are released on 3 conveyors, arranged in unit loads and wrapped in so-called columns of pallets. Conveyors hold a limited number of columns (precisely, 10, 14 and 8 columns, respectively) and need to be emptied as soon as possible not to block subsequent releases. The storing list contains the details of each release. In particular, separately per the product type, it reports the release time, the amount of columns released per the product type, the conveyor of release as well as the stacks assigned to the product type for storing, and the order of precedence according to which they have to be filled up. The decisions on assignment and sequencing of stacks per product type are taken according to Lanza et al. (2022a, 2022b).

Orders are shipped during the first and the second shift only, while the third shift is dedicated to filling the collection area as much as possible to quickly begin the truck loading operations the next morning. The shipping list of the day is generally known a day in advance and describes the composition of each order, that is, types of product and amount of columns requested, as well as the leaving time of the associated truck. Items are required to be retrieved from stacks following the given order of precedence per product type and to be moved to the collection area before a given due date, not to generate truck loading delay.

The fleet of the company is composed of 5 LGV autonomous shuttles and 7 human-driven forklifts, some of which are electric and equipped with a lithium-ion battery. By considering the general problem description in Sect. 3, they correspond to vehicles of type F1 and type F2, respectively. In the following, we will refer to the above-mentioned types of vehicles as LGV and FKL for short. Both types of vehicles may transport 2 columns at most at the same time, starting from and ending at their parking slot. The parking slot of the FKL vehicles also hosts the charging station.

The structure of the warehouse and the internal layout of each department are analogous to the ones shown in Fig. 1 in Sect. 3.

6.2 Plan of the experiments

We report the results of three types of experiments, as detailed below.

• 1. In the first type of experiments, we evaluate the matheuristic approaches presented in Sect. 5 on a set of small-medium size instances, for which the state-of-the-art optimization solver CPLEX was able to determine optimal solutions or lower bounds (Sect. 6.4.1). In fact, CPLEX did not succeed in solving the model in Sect. 4.3 on real-size instances, due to their dimension.

• 2. In the second type of experiments, we provide some technical and managerial insights by using the set of small-medium size instances considered in the first type of experiments (Sect. 6.4.2).

• 3. In the third type of experiments, we investigate the efficiency and the efficacy of algorithm ATD, which proved to be the most promising matheuristic according to the first type of analysis, on a wide pool of real instances related to the addressed case study (Sect. 6.5).

Both the small-medium size instances and the real instances are described in Sect. 6.3. The formulation as well as the matheuristic approaches have been implemented using the OPL language and solved via CPLEX 12.9 solver (IBM ILOG, 2016). All the experiments have been conducted on an Intel Xeon 5120 computer with 2.20 GHz and 32 GB of RAM.

6.3 The instances

The reference company provided us with the following information for a pool of selected 8-hour shifts:

• (i) the warehouse configuration at the beginning of the shift, i.e., product types and corresponding number of columns inside the warehouse;

• (ii) the storing list of the shift;

• (iii) the shipping list of the shift and of the next 3 shifts.

Some data needed to be integrated, while others were randomly generated since not provided by the company. Specifically, the positions of the columns in the warehouse at the beginning of a shift are randomly generated by respecting some agreed industrial practice or insights provided by the company, to have a realistic configuration. Also, we randomly generated the retrieval order of precedence per product type for the occupied stacks. Regarding the storing list, the set of stacks assigned to each product type and the corresponding filling order of precedence were obtained by the method in Lanza et al. (2022a). While the number of columns released during the shift is provided, release time instants have been randomly generated. Concerning the shipping list, the composition of each order is given. However, the truck leaving times have been randomly generated by considering that the majority of the orders are shipped in the morning.

In the first type of experiments, we tested 15 instances of small-medium size, which have been generated starting from the above-described real data set. To be able to find optima with CPLEX, in such instances, the duration of a shift is shortened from 8 to 4 h, and the number of product types and columns to move is reduced proportionally to the shortening of the planning horizon. On the other hand, the number of vehicles used to perform the operations, the number of conveyors and the number of collectors are unchanged.

The features of such 15 instances are described in Table 2, which reports the number of storage locations, SL, together with the number of product types in ${\mathcal {K}}_{in}$ and in ${\mathcal {K}}_{out}$ (columns $K_{in}$ and $K_{out}$, respectively) and the corresponding number of items to move (columns $C_{in}$ and $C_{out}$, respectively). This information has been also summarized in Table 2a. The instances are partitioned into 3 classes. In the first class, $C_{in}$ and $C_{out}$ are about the same (Table 2b), in the second class, $C_{in}$ is higher than $C_{out}$ (Table 2c), and in the third class, $C_{in}$ is lower than $C_{out}$ (Table 2d).

Table 2 Features of small-medium size instances

For the third type of experiments, instead, real shifts have been considered, thus generating corresponding real instances. In particular, 15 instances corresponding to 5 ordinary days are first considered (i.e., 5 days times 3 daily shifts). Moreover, the other 15 instances, corresponding to 5 very busy days with respect to both production and shipments, are taken into account. In such 5 very busy days both production and shipments are higher at about 25% with respect to a normal day, and about 500 more movements are required for storing or retrieving items per shift. In these experiments, the 15 instances of each group are solved in cascade, from the instance corresponding to the first shift of the first day, to the instance corresponding to the last shift of the fifth day. The features of the real instances are reported in Table 3 separately for the first and the second group of days, by reporting the same kind of information summarized in Table 2a.

Table 3 Features of real instances

Regarding the parameters related to the battery management, $B^-$, $B^+$ and $\Theta$ have been set to $30\%$, $80\%$ and $35\%$ of the total capacity of the battery, respectively. In particular, $B^-$ and $B^+$ are set following the features described in the user manual of the specific type of electric FKL, denoted as eFKL in the following, used by the company. On the other hand, the parameters $e^{ij}$, $e^{ijk}$ and $e^k$ have been calculated following the procedure described in Sect. 4.2 and using the data in Table 4. Specifically, Table 4 reports some real data related to the type of forklift used by the reference company in the considered warehouse, i.e., a STILL RX-20 electric forklift truck whose technical data sheet is also available on the web. For those data not available, we performed a reasonable approximation by taking into account the physical characteristics of the forklift. Finally, the value of the parameters $c_{max}$, $\psi ^{v}_0$ and $\Delta$ does depend on the kind of analysis performed, and therefore it will be specified in the following sections.

Table 4 Energy consumption model parameters

6.4 Computational results on small-medium instances

6.4.1 Matheuristics evaluation

As outlined, the first type of experiments aims to evaluate the performance of the matheuristics described in Sect. 5 on the small-medium size instances of our data set, for which an optimal value, or a lower bound, can be obtained via CPLEX.

The proposed formulation relies on parameters $\psi$, $\xi$ and $\zeta$, which are related to the soft optimization criteria. In particular, increasing values of $\psi$ would give priority to emptying conveyors and moving columns as soon as they are released from the production area, while increasing values of $\xi$ would give priority to the anticipation moves toward the collection area. Parameter $\zeta$, instead, controls the charging process, trying to limit useless (and costly) charges. An analysis of the first two parameters can be found in Lanza et al. (2022d). The results reported in the following are obtained by setting the value of both parameters to 10 since this combination proved to be very effective. Regarding $\zeta$, after some preliminary tests, it has been set to 1, so as not to discourage the use of the EVs and, at the same time, limit the useless recharge of the batteries.

We solved the 15 small-medium size instances described in Sect. 6.3 with the 3 matheuristics. Each instance has been solved 3 times, by varying the number of the EVs. Specifically, 1, 2 and 3 eFKLs have been considered (recall that the total number of FKLs is 7). After some preliminary tests, we considered 3 different settings for the initial level of charge $\psi ^v_0$ of the eFKLs depending on the number of EVs. Specifically, if one eFKL is used, its initial level of charge is set to half of the range $[B^-,B^+]$; if 2 eFKLs are used, the range $[B^-,B^+]$ is split into 2 parts of equal length, and the initial level of charge of one vehicle is set to the half of the first range, while the initial level of charge of the other vehicle is set to the half of the second range; finally, if 3 eFKLs are used, then $[B^-,B^+]$ is split into 3 parts of equal length, and the initial level of charge of the vehicles is set to the half of the first, of the second and of the third range, respectively. Finally, the tests have been performed by setting $c_{max} = 1$, i.e., only 1 eFKL can be charged at a time, and $\Delta = 0$. Note, however, that parameter $\Theta$ still regulates the level of charge of the eFKLs at the end of the time horizon, by guaranteeing a minimum level of charge to each of them.

Now we describe the parameter setting of the three matheuristics. For TD, we split the 4 h time horizon into 4 subshifts, thus obtaining subshifts of about 60 min. As reported in Lanza et al. (2022d), longer subshifts may lead to hardly solvable subproblems, while shorter subshifts seem to negatively affect the quality of the solutions obtained. The resolution of each subproblem has been performed via CPLEX by stopping the execution as soon as an optimality gap less than $3\%$ or a time limit of 15 min was reached. Most subproblems, however, were solved to optimality. For AD, both phases are solved via CPLEX and phase 1 (i.e., the resolution of $P_{K_{in}}$) is stopped as soon as the optimality gap is less than 1%, while phase 2 (i.e., the resolution of $P_{K_{out}}$) is stopped as soon as the optimality gap is less than 3%. Alternatively, for both phases, a time limit of 30 minutes is imposed. A finer optimality gap for phase 1 has been chosen for two reasons. Firstly, the resolution of $P_{K_{in}}$ is more critical and solutions of poor quality may jeopardize the anticipation moves, which are planned in phase 2; moreover, the resolution of $P_{K_{in}}$ is faster and good optimality gaps are more likely to be obtained within the time limit. Finally, for ATD we split the time horizon into 4 subshifts. In each subshift, we stop phase 1 (i.e., the resolution of $P_{K_{in}}$) as soon as the optimality gap is less than 1%, and phase 2 (i.e., the resolution of $P_{K_{out}}$) as soon as the optimality gap is less than 3%, as for AD. Alternatively, for both phases, when a time limit of 7 minutes is reached.

For each instance, Tables 5, 6 and 7 report the solving times and the optimality gaps related to the 3 matheuristics, for each of the 3 fleet compositions mentioned before. For each matheuristic, the solving time (expressed in seconds) is the sum of the times needed to solve the related subproblems. On the other hand, the optimality gaps are calculated with respect to the optimal value found by CPLEX by solving GSRP with only 1 eFKL. This is the only scenario of GSRP that CPLEX was able to solve to optimality and the corresponding optimal values thus represent lower bounds in the case of 2 and 3 eFKLs. The times required to compute such lower bounds are reported in column LB. Finally, the last 4 rows of each table report average values: in the first 3 rows, the averages refer to the three classes in which the 15 instances are partitioned, namely $C_{in} \simeq C_{out}$, $C_{in} > C_{out}$, and $C_{in} < C_{out}$, while the last row reports the average over the 15 instances. If an instance could not be solved within the time limit, then the symbol ‘-’ appears. In this case, the averages are calculated by considering only the solved instances.

Table 5 Performance of the matheuristics for 1 eFKL

Table 6 Performance of the matheuristics for 2 eFKLs

Table 7 Performance of the matheuristics for 3 eFKLs

In the case of 1 eFKL, CPLEX can find optimal solutions to all 15 instances, with an average computational time of half an hour (note, however, that some instances required more than 1 h to be solved). The third class, in which fewer columns of type $K_{in}$ need to be moved than columns of type $K_{out}$, i.e., $C_{in} < C_{out}$, seems to be the most tricky to address. All the three matheuristics, i.e., TD, AD and ATD, show a very good performance, being able to solve the instances in a much shorter time than CPLEX (see column LB), as expected, with an average optimality gap of $5.70\%$, $8.48\%$, and $2.61\%$, respectively. Note however that only AD and ATD found one optimal solution. TD and AD seem to have opposite behaviour, the first one outperforming the second one with respect to the average optimality gap, while the second one outperforming the first one with respect to average computational time. In particular, TD shows a better average optimality gap for the class of instances where more columns of type $K_{in}$ need to be moved than columns of type $K_{out}$ (i.e., $C_{in} > C_{out}$). As opposed, AD seems to have a more stable behaviour among the three classes of instances. Finally, as mentioned, AD can find an optimal solution for instance 2. ATD improves the good features of both TD and AD. It outperforms TD and AD in both average computational time and average optimality gap, being still able to find an optimal solution for instance 2. Better results in terms of average optimality gap are obtained for the class $C_{in} > C_{out}$, while for the other two classes, better results are obtained on average for the computational time. Compared to CPLEX, the average reduction of computational time is $90\%$.

When 2 eFKLs are used, AD can solve only 9 out of 15 instances, being still able to find an optimal solution to instance 2. On the other hand, TD can find a feasible solution for all 15 instances. Both computational time and optimality gap of TD are outperformed on average by ATD, which can find better solutions, in terms of average optimality gap, almost one-third of the time on average. ATD shows more stable computational times on average for the three classes of instances, still showing slightly better performance in terms of average optimality gap for the class $C_{in} > C_{out}$. It is also able to find an optimal solution for instance 2.

Finally, with 3 eFKLs the instances become much more tricky to address. AD is not able to find any feasible solution within the time limit imposed, while TD solves 10 instances out of 15. On the other hand, ATD still shows very good behaviour, being able to find solutions to all the instances with an average gap of about 5% in a short time. Interestingly, differently than in the scenarios with 1 or 2 eFKLs, the class $C_{in} > C_{out}$ seems now the most tricky to address.

To investigate the performance of TD, AD and ATD more accurately, we also report some performance profiles. The performance profile for a solver is the cumulative distribution function of a performance measure (see Dolan and Moré 2002). The comparison of performance profiles of different solvers may provide useful information about the relative performance of one solver against the others, often hidden when only comparing average results. Figure 3 reports the performance profiles of TD, AD and ATD related to the solving time (Fig. 3a) and the optimality gap (Fig. 3b), only considering the case of 1 eFKL, for which all the approaches can solve all the 15 instances. According to such profiles, the dominance of ATD over TD and AD emerges in terms of the optimality gap. Interestingly, AD solves all the instances in less than 500 s, generally outperforming both TD and ATD in terms of computational time.

Fig. 3

Performance profiles of TD, AD and ATD with respect to solving time and optimality gap in the case of 1 eFKL

To better analyse the results in Tables 5, 6 and 7, and to provide some insights on the quality of the computed solutions, we report in Tables 8, 9 and 10 some aggregated features of the computed solutions in terms of crucial performance indicators suggested by the reference company, by considering the 3 matheuristics and the 3 fleet compositions. Indeed, for each fleet composition, only the approaches to solving all 15 instances are considered. Specifically, the primary goal is analysed in terms of the average time, in minutes, travelled by an LGV and by an FKL over the 15 instances (LGV Avg. Travel Time and FKL Avg. Travel Time, respectively). The results are also reported separately for conventional and electric FKLs (Conventional FKL Avg. Travel Time and Electric FKL Avg. Travel Time, respectively), calculating the averages over the corresponding number of vehicles. Moreover, we report the average charging time, in minutes, of an eFKL (Electric FKL Avg. Charging Time). The other secondary goals, i.e., emptying conveyors and anticipation moves, are evaluated by considering the average time, in minutes, an incoming item idles on an input point before being moved to an available collector (Input point Avg. Idle Time per item), and the average percentage of columns, requested in the next shift, moved to the collection area at the end of the planning horizon (Saturation of collection area). In the case of 1 eFKL, the average time travelled by an LGV is pretty similar in the solutions obtained by the 3 matheuristics, with ATD slightly outperforming such an indicator with respect to AD and TD. On the other hand, the average time travelled by an FKL is the highest for ATD with respect to the other two approaches. TD shows a longer average travel time for the eFKL, and consequently longer average charging time. The lowest average travel time for the eFKL is achieved by AD for which, accordingly, also the lowest average charging time is obtained. Regarding the soft objectives, very good results are reported for all the approaches. The lowest average idle time on conveyors for incoming items is obtained by ATD. This may be explained by considering the structure of the algorithm, which gives priority to storing operations (being $P_{K_{in}}$ solved first) at the expense of the picking ones (being $P_{K_{out}}$ dependent on the solution of $P_{K_{in}}$). Additionally, the subproblems to handle are easier to solve with respect to AD (as already observed in the previous analysis, see Table 5), and therefore better performance can be observed with respect to AD, whose subproblems $P_{K_{in}}$ are still tricky to address within the time limit. Nevertheless, AD has the best results for the anticipation moves, being able to move to the collection area the 98.67% of the columns requested in the next shift. This may be explained by observing that AD has a complete view of the planning horizon, as opposed to TD and ATD.

Table 8 Features of solutions of TD, AD and ATD for 1 eFKL

Table 9 Features of solutions of TD and ATD for 2 eFKLs

Table 10 Features of solutions of ATD for 3 eFKLs

With 2 eFKLs, pretty similar results may be observed for TD and ATD in terms of the average travel time of LGV and FKL. ATD shows a shorter average travel time for the eFKLs with, however, higher charging time on average. Moreover, ATD outperforms TD for the other two secondary goals, slightly worsening the saturation of the collection area with respect to the case of 1 eFKL.

Finally, with 3 eFKLs the average travel time of LGV and FKL slightly increases with respect to the scenarios with 1 or 2 eFKLs. The increased average travel time of the FKLs, however, is mainly due to the intensification of activities assigned to the eFKLs, thus promoting sustainability in operations. In fact, the average time travelled by a traditional FKL decreases when the number of eFKLs increases. Consequently, longer recharges are required. The more frequent recharges slightly affect the average travel time of an eFKL, which is higher compared to the previous two cases since the eFKL needs to reach the charging station more often. In spite of this, the average idle time of incoming items on conveyors is the same as for 1 and 2 eFKLs, while the percentage saturation of the collection area at the end of the time horizon is only negligibly worse.

According to the results presented above, the matheuristic ATD appears to be very promising in all the tested scenarios, both in terms of solving time and optimality gap, usually outperforming TD and AD, and also considering relevant quality indicators like the ones previously commented.

6.4.2 Technical and managerial insights

In this set of experiments, focusing on the green aspects of the problem, we performed an analysis aimed at deriving both some technical insights on what makes GSRP hard to solve, and also some managerial insights related to the quality indicators reported in Sect. 6.4.1.

Regarding the technical analysis, the following aspects are considered: the number of eFKLs (parameter $|{\mathcal {V}}^{E}|$), the minimum percentage that the total charge of all batteries, at the end of the time horizon, is imposed to have with respect to the initial total charge of the batteries (parameter $\Delta$), and the maximum number of EVs that can be charged simultaneously at the charging station (parameter $c_{max}$). As indicated before, we considered 3 fleet compositions, namely with 1, 2 or 3 eFKLs. For each of the 3 fleet compositions, we investigated 4 values of parameter $\Delta$, namely $0.25, \, 0.5, \, 0.75, \, 1$. Additionally, for each fleet composition and each value of $\Delta$, we considered both the case in which only 1 eFKL at a time may be charged, and the case where all the eFKLs composing the fleet may be charged simultaneously. The analysis has been performed by considering the matheuristic ATD, which proved to be very efficient in solving GSRP on instances of small to medium size (see Sect. 6.4.1). For each of the above-mentioned settings, in Table 11 we report the number of instances (out of 15) that ATD was able to solve within the time limit (i.e., 7 min for each phase), and the average time, in seconds, required by ATD to solve such instances.

Table 11 Some technical insights by varying battery management features

Let us first consider the impact of parameter $\Delta$. By increasing $\Delta$, the instances become trickier to handle by ATD, for all the 3 fleet compositions. In particular, the number of instances solved by ATD decreases for each setting of $|{\mathcal {V}}^{E}|$ and $c_{max}$. For $|{\mathcal {V}}^{E}|=1$, ATD is able to solve all the instances for $\Delta = 25\%$ and $\Delta = 50\%$. Nevertheless, its solving time increases by about $200\%$. The lower number of instances are solved for $|{\mathcal {V}}^{E}|=3$ and $c_{max}=1$, which is the hardest setting. In particular, no instance is solved for $\Delta = 100\%$.

Let us consider now the impact of the number of the eFKLs when only one vehicle at a time may be charged at the charging station (i.e., $|{\mathcal {V}}^{E}|=1,2,3$ and $c_{max}=1$ in Table 11). An increased number of eFKLs makes the instances more difficult to solve. In fact, more EVs are present in the fleet, fewer instances are solved by ATD. For $\Delta = 25\%$, such a number decreases from 15 instances, in the case of 1 eFKL, to 11 instances, in the case of 3 eFKLs. A more significant reduction of solved instances can be observed for $\Delta = 50\%$. Finally, when $\Delta =100\%$ the reduction of solved instances is remarkable, passing from 11 in case of $|{\mathcal {V}}^{E}|=1$, to 2 and 0 for $|{\mathcal {V}}^{E}|=2$ and $|{\mathcal {V}}^{E}|=3$, respectively.

Finally, some benefits may be observed in terms of the number of solved instances when the simultaneous charge of more eFKLs is allowed. Let us compare $|{\mathcal {V}}^{E}|=2$ and $c_{max}=1,2$, with $|{\mathcal {V}}^{E}|=3$ and $c_{max}=1,3$. For $|{\mathcal {V}}^{E}|=2$ and $\Delta = 25\%$, a simultaneous recharge allows one to solve all the 15 instances, compared to the 12 solved instances in case of $c_{max}=1$. Even, in the worst setting, i.e., $|{\mathcal {V}}^{E}|=3$ and $\Delta = 100\%$, the simultaneous recharge allows one to solve 3 instances instead of none of them, as happens when $c_{max}=1$.

By summarizing, it is possible to outline the following technical insights:

• increasing the value of parameter $\Delta$ makes the instances harder to solve; in particular, $\Delta = 25\%$ seems to be a good option;

• increasing the number of eFKLs makes instances more difficult to solve as well;

• with simultaneous recharge, i.e., $c_{max} \ge 2$, ATD is able to solve more instances with respect to the scenario where only 1 eFKL at a time may recharge.

Lastly, by considering the scenario $\Delta = 25\%$, for which more instances are solved according to  Table 11, we analyse the quality performance indicators suggested by the reference company (the same discussed in Sect. 6.4.1) for the subset of instances solved by ATD for all the combinations $|{\mathcal {V}}^{E}|= 1, \, 2, \, 3$ and $c_{max}=1, 2, 3$. The number of instances that are solved by ATD for all the above-mentioned combinations of $|{\mathcal {V}}^{E}|$ and $c_{max}$ is 10, thus the averages reported in Table 12 are calculated over 10 instances. We can observe that the average time travelled by an LGV is pretty similar in all the considered scenarios. As opposed, the average time travelled by an FKL is more floating, and it depends on the considered green features. Focusing on such average travelled times separately for eFKL and for conventional FKL, the highest average travel time for eFKL is observed for $|{\mathcal {V}}^{E}|= 2$ and $c_{max}=1$, and for $|{\mathcal {V}}^{E}|= 3$ and $c_{max}=3$, with a correspondent reduction of the time travelled by a conventional FKL. The selection of such parameter settings seems thus to be able to enhance the green management of the warehouse, therefore improving sustainability and healthiness for workers, e.g., in terms of reduced noise and better local air quality. Regarding the average charging time, the highest values are observed when 3 eFKLs are present in the fleet, and this is quite expected due to the larger amount of operations performed by the EVs in such scenarios. Such a larger charging time may determine increased costs for recharging the EVs, and consequently an increased management cost for the warehouse. Notice, however, that thanks to the limitation of the unnecessary charging periods modelled through the last term of the objective function (5), such costs are indeed the ones which are essential to achieve the enhanced level of sustainability within the warehouse. Finally, regarding the other soft objectives pursued, i.e., average idle time of incoming items on conveyors and saturation of the collection area at the end of the time horizon, very good results are overall achieved. Notice, in particular, that anticipation moves appear to be particularly effective when more EVs are used.

The performed tests, although very preliminary, seem thus to suggest that using more EVs in the fleet, and calibrating the maximum number of vehicles to be recharged simultaneously as a function of their number, may enhance the sustainability and, in general, the green management of the warehouse, still maintaining very good values of relevant quality performance indicators, like the idle time of the incoming items on the conveyors, and the saturation of the collection area at the end of the time horizon, thanks to the performed anticipation moves.

Table 12 Some managerial insights by varying battery management features

6.5 Computational results on real instances

We tested the efficiency of ATD on a real testbed composed of 15 instances corresponding to 5 ordinary days of work (remember that a day comprises 3 shifts, leading to 3 instances), and 15 instances corresponding to 5 extremely busy days of work, where both storing and retrieving operations are particularly crucial to manage, just before a peak period of requests. We always considered 1 eFKL. Moreover, to adhere as much as possible to the realistic scenario where operations are performed 24 h, we set $\Delta = 1$, i.e., the total level of charge of the eFKLs at the end of each shift must be at least equal to the initial level of charge.

ATD was able to determine a solution to all the shifts composing the 10 days, by still maintaining good performance for the company indicators. Table 13 reports the same indicators of Table 5 except for the gap since the real instances could not be solved by CPLEX. In the first column (Ord. days), average results related to the ordinary days are reported, while in the second one (Extr. busy days) the average results concern the extremely busy days. Finally, in the third column, we report the percentage difference between the first and the second column.

The average computational time required by ATD to solve busy days is much higher than the time to solve ordinary days (+87%). However, it is still lower than the time required by the company to solve an entire shift, which is four hours. The increase in solving time for busy days is especially due to the greater difficulty of ATD in solving the $P^{\lambda }_{K_{in}}$ subproblems. The increased number of movements on the selected busy days causes an unavoidable increase in travel times for both LGVs and FKLs when compared to the ordinary days (+27% and +19%, respectively). The increase in travel time seems to equally affect both types of FKLs, i.e., conventional and electric ones (in both cases the increase in travel time is $+19\%$). The longer travel time for the eFKL during busy days causes a consequent greater need for recharging ($+15\%$) compared to ordinary days. A slightly longer idle time for items on conveyors is observed on busy days. Finally, the saturation of the collection area does not seem to be penalized on busy days, being this indicator almost the same in both cases.

Figures 4 and 5 report the trend of two crucial features of the problem, for both the 5 ordinary and the 5 busy days. Specifically, Figs. 4a and 5a report the level of charge of the battery of the eFKL, while Figs. 4b and 5b report the saturation trend of the collection area. Note that in both datasets there are shifts where the collection area is empty and not replenished (see the last shift of day 3 in Fig. 4b, and the last shift of day 5 in Fig. 5b). The motivation is that the shipment list of the subsequent day is not yet available, and anticipation moves cannot be performed. Disregarding such days in the computation of the ‘Saturation of collection area’ of Table 13, the values become $84\%$ and $83\%$, thus testifying the very good exploitation of the collection area. Figures 4a and 5a, instead, show the typical trend for the level of charge of a battery allowing partial recharging, which can be continuously recharged during operational downtime. In conclusion, ATD appears to be a valuable approach for solving GSRP also in real scenarios.

Table 13 Features of solutions of ATD in terms of crucial performance indicators, for ordinary and busy days

Fig. 4

Trends for ordinary days solved with ATD

Fig. 5

Trends for busy days solved with ATD

7 Conclusions

We considered the green sequencing and routing problem (GSRP), which is a generalization of the Sequence and Routing problem consisting in defining the best sequence of locations to visit within a warehouse, for storing and/or retrieving items in a given time horizon, with a fleet composed of both electric and conventional vehicles. Firstly, we generalized a mathematical model already proposed by the authors with additional battery management features. Then, we proposed two matheuristic algorithms for GSRP: the first one is based on activity decomposition, while the second one combines activity decomposition with time horizon decomposition, the latter being the core of a different matheuristic previously developed by the authors. The matheuristics have been evaluated on a set of small-medium instances showing very good results both in terms of algorithmic performance and solution quality and proving the superiority of the combined activity and time horizon decomposition approach. The combined approach has been tested also on two sets of real instances, corresponding to ordinary days and extremely busy days of work, respectively, showing very good results also in such real scenarios. Some technical insights on what makes GSRP hard to solve, and some managerial insights focusing on the green aspects of the problem, have been also provided.

Future research will investigate alternative battery management policies such as the swapping of batteries, instead of partial recharging as analysed in this paper. Battery swapping introduces a range of different challenges in comparison to partial recharging for the operation of EVs in warehouses. These challenges include not only defining a suitable battery inventory but also meticulously scheduling recharging procedures for depleted batteries to prevent disruptions to EV operations. Additionally, synchronizing battery swapping activities with the schedules of single EVs adds another layer of complexity. To address this challenging problem, we intend to develop a Mixed-Integer Linear Programming formulation aimed at addressing the distinctive aspects of this problem. Furthermore, recognizing the crucial role that product placements within the warehouse play in influencing vehicle routes, we plan to extend this study by integrating storage location assignment into the GSRP framework. Such integration should yield even more optimized solutions in terms of energy consumption within the warehouse, promoting environmentally sustainable operations. Nevertheless, given the inherent complexity of this integrated problem (both the storage location assignment problem and the GSRP fall within the NP-hard class), we plan to create specific heuristic or matheuristic algorithms to tackle realistic instances.

Appendix

The notation used in the mathematical formulation is summarized in Table 14. We denote by ${\mathcal {N}}^+(i)$ and ${\mathcal {N}}^-(i)$ the sets of nodes linked to $i \in \mathcal {N}$ via an exiting and an entering arc, respectively, that is

$$\begin{aligned} {\mathcal {N}}^+(i) = \left\{ j \in {\mathcal {N}}:\ \exists \ (i,j) \in {\mathcal {A}}\right\} , \qquad {\mathcal {N}}^-(i) = \left\{ j \in {\mathcal {N}}:\ \exists \ (j,i) \in {\mathcal {A}}\right\} . \end{aligned}$$

Table 14 Sets, parameters and variables used in the model

Additional constraints to formulation (5)–(16):

$$\begin{aligned}&\begin{aligned} \sum _{(j,t') \in {\mathcal {N}}^{+}(i,t)} x^v_{(i,t)(j,t')} - \sum _{(j,t') \in {\mathcal {N}}^{-}(i,t)} x^v_{(j,t')(i,t)} = {\left\{ \begin{array}{ll} 1 &{}\text{ if } (i,t) = (\omega ^{1},\,0), \\ -1 &{} \text{ if } (i,t) = (\omega ^{1},\,T), \\ 0 &{}\text{ otherwise }, \end{array}\right. } \\ \forall \ (i,t) \in {\mathcal {N}}_{F1}, \ \forall \ v \in {\mathcal {V}}^{1}, \end{aligned}\end{aligned}$$

(18)

$$\begin{aligned}&\begin{aligned} \sum _{(j,t') \in {\mathcal {N}}^{+}(i,t)} x^v_{(i,t)(j,t')}-\sum _{(j,t') \in {\mathcal {N}}^{-}(i,t)} x^v_{(j,t')(i,t)} = {\left\{ \begin{array}{ll} 1 &{}\text{ if } (i,t) = (\omega ^{2},\,0), \\ -1 &{} \text{ if } (i,t) = (\omega ^{2},\,T), \\ 0 &{}\text{ otherwise }, \end{array}\right. } \\ \forall \ (i,t) \in {\mathcal {N}}_{F2}, \forall \ v \in {\mathcal {V}}^{2}, \end{aligned} \end{aligned}$$

(19)

$$\begin{aligned}&\sum _{v \in {\mathcal {V}}^{1}} x^v_{(i,t)(j,t')} \le 1&\forall \ ((i,t),(j,t')) \in {\mathcal {A}}_{F1}:\ i, j \ne \omega ^{1}, \end{aligned}$$

(20)

$$\begin{aligned}&\sum _{v \in {\mathcal {V}}^{2}} x^v_{(i,t)(j,t')} \le 1&\forall \ ((i,t),(j,t')) \in {\mathcal {A}}_{F2}:\ i,j \ne \omega ^{2}, c, \end{aligned}$$

(21)

$$\begin{aligned}&\begin{aligned}&\sum _{(j,t') \in {\mathcal {N}}^{+}(i,t)} y^k_{(i,t)(j,t')}-\sum _{(j,t') \in {\mathcal {N}}^{-}(i,t)} y^k_{(j,t')(i,t)} \\&\qquad = {\left\{ \begin{array}{ll} d_{in}^k(i,t) + u^k_{i} &{}\text{ if } i \in {\mathcal {R}}, \, t = 0, \\ u^k_{i} &{}\text{ if } i \in {\mathcal {B}}, \, t = 0, \\ d_{in}^k(i,t) &{} \text{ if } i \in {\mathcal {R}}, \, t = 1, \ldots , T, \\ 0 &{} \text{ if } i \in {\mathcal {B}}, t = 1, \ldots , T, \\ 0 &{}\text{ if } i \in {\mathcal {S}}^{k}_{out} \cup {\mathcal {S}}^{k'}, \, t = 0, \ldots , T, \end{array}\right. } \\&\forall \, k \in {\mathcal {K}}_{in},\, \forall \, k' \in {\mathcal {K}}:\ k' \ne k, \\&\forall \, (i,t) \in {\mathcal {N}}_{in}:\ i \in {\mathcal {R}} \cup {\mathcal {B}} \cup {\mathcal {S}}^{k}_{out} \cup {\mathcal {S}}^{k'}, \end{aligned} \end{aligned}$$

(22)

$$\begin{aligned}&\begin{aligned} \sum _{(j,t') \in {\mathcal {N}}^{+}(i,t)} y^k_{(i,t)(j,t')} -\sum _{(j,t') \in {\mathcal {N}}^{-}(i,t)} y^k_{(j,t')(i,t)} = {\left\{ \begin{array}{ll} u^k_{i} &{} \text{ if } i \in {\mathcal {B}}, t = 0, \\ 0 &{} \text{ if } i \in {\mathcal {B}}, t \ge 1, \\ 0 &{}\text{ if } i \in {\mathcal {S}}^{k}_{in} \cup {\mathcal {S}}^{k'} \\ &{} \text {and }t \ge 0, \end{array}\right. } \\ \forall \ k \in {\mathcal {K}}_{out}, \ \forall \ k' \in {\mathcal {K}}:\ k \ne k', \\ \forall \ (i,t) \in {\mathcal {N}}_{out}:\ i \in {\mathcal {B}} \cup {\mathcal {S}}^{k'} \cup {\mathcal {S}}^{k}_{in}, \end{aligned} \end{aligned}$$

(23)

$$\begin{aligned}&\begin{aligned} \sum _{(j,t') \in {\mathcal {N}}^{-}(\pi ,t)} y^k_{(j,t')(\pi ,t)} - \sum _{(\pi ,t') \in {\mathcal {N}}^{+}(\pi ,t)} y^k_{(\pi ,t)(\pi ,t')} = d^k_{out}(\pi ,t) \\ \forall \ k \in {\mathcal {K}}_{out}, \ \forall \ t \ge 1, \end{aligned} \end{aligned}$$

(24)

$$\begin{aligned}&y^k_{(\pi ,0)(\pi ,1)} = u^k_{\pi } \qquad \forall \ k \in {\mathcal {K}}_{out}, \end{aligned}$$

(25)

$$\begin{aligned}&\sum _{\begin{array}{c} k \in {\mathcal {K}}_{in}: \\ ((i,t),(j,t')) \in {\mathcal {A}}_{in} \end{array}} y^k_{(i,t)(j,t')} \le c_{F1} \sum _{v \in {\mathcal {V}}^{1}} x^v_{(i,t)(j,t')} \qquad \forall \ ((i,t),(j,t')) \in {\mathcal {A}}_{F1}, \end{aligned}$$

(26)

$$\begin{aligned}&\sum _{\begin{array}{c} k \in {\mathcal {K}}_{in}: \\ ((i,t),(j,t')) \in {\mathcal {A}}_{in} \end{array}} y^k_{(i,t)(j,t')} + \sum _{\begin{array}{c} k \in {\mathcal {K}}_{out}: \\ ((i,t),(j,t')) \in {\mathcal {A}}_{out} \end{array}} y^k_{(i,t)(j,t')} \le c_{F2} \sum _{v \in {\mathcal {V}}^{2}} x^v_{(i,t)(j,t')}\nonumber \\&\quad \qquad \qquad \qquad \quad \qquad \qquad \qquad \quad \qquad \qquad \qquad \quad \qquad \qquad \forall \ ((i,t),(j,t')) \in {\mathcal {A}}_{F2}, \end{aligned}$$

(27)

$$\begin{aligned} \sum _{k \in {\mathcal {K}}_{in}} d_{in}^k(r,t) + \sum _{k \in {\mathcal {K}}_{in}} y^k_{(r,t-1)(r,t)}&\le {\left\{ \begin{array}{ll} \displaystyle c_{r} - \sum _{k \in {\mathcal {K}}_{in}} u^k_{r} &{}\text {if }t = 1, \\ c_{r} &{}\text{ if } t > 1, \end{array}\right. }\nonumber \\&\quad \quad \qquad \forall \ r \in {\mathcal {R}},\ \forall \ t \ge 1, \end{aligned}$$

(28)

$$\begin{aligned} \begin{aligned} \sum _{(j,t') \in {\mathcal {N}}^{-}(b,t)} \quad \sum _{k \in {\mathcal {K}}} y^k_{(j,t')(b,t)} \le {\left\{ \begin{array}{ll} \displaystyle c_b - \sum _{k \in {\mathcal {K}}} u^k_{b} &{}\text {if }t = 1,\\ c_{b} &{}\text{ if } t > 1, \end{array}\right. } \\ \forall \ b \in {\mathcal {B}}, \ \forall \ t \ge 1, \end{aligned} \end{aligned}$$

(29)

$$\begin{aligned}&\sum _{(j,t') \in {\mathcal {N}}^{-}(\pi ,t)} \quad \sum _{k \in {\mathcal {K}}_{out}} y^k_{(j,t')(\pi ,t)} \le {\left\{ \begin{array}{ll} \displaystyle c_{\pi } - \sum _{k \in {\mathcal {K}}_{out}} u^k_{\pi } &{}\text{ if } t = 1, \\ c_{\pi } &{}\text{ if } t > 1, \end{array}\right. } \quad \forall \ t \ge 1, \end{aligned}$$

(30)

$$\begin{aligned}&\begin{aligned} \sum _{\tilde{t}=0}^{t} \left[ \sum _{(j,t') \in {\mathcal {N}}^{-}(i,\tilde{t})} y^k_{(j,t')(i,\tilde{t})} - \sum _{(j,t') \in {\mathcal {N}}^{+}(i,\tilde{t})} y^k_{(i,\tilde{t})(j,t')} \right] \le c_{i} \\ \forall \ k \in {\mathcal {K}}_{in}, \ \forall \ t \ge 1, \ \forall \ (i,t) \in \mathcal {N}_{in}:\ i \in {\mathcal {S}}^{k}_{in}, \end{aligned} \end{aligned}$$

(31)

$$\begin{aligned}&\begin{aligned} \sum _{\tilde{t}=0}^{t} \left[ \sum _{(j,t') \in {\mathcal {N}}^{+}(i,\tilde{t})} y^k_{(i,\tilde{t})(j,t')} -\sum _{(j,t') \in {\mathcal {N}}^{-}(i,\tilde{t})} y^k_{(j,t')(i,\tilde{t})} \right] \le c_{i} \\ \forall \ k \in {\mathcal {K}}_{out}, \ \forall \ t \ge 1, \ \forall \ (i,t) \in \mathcal {N}_{out}:\ i \in {\mathcal {S}}^{k}_{out}, \end{aligned} \end{aligned}$$

(32)

$$\begin{aligned}&\sigma ^t_{s^k} = \sum _{\tilde{t}=0}^t \left[ \sum _{(j,t') \in {\mathcal {N}}^{-}(s^k,\tilde{t})} y^k_{(j,t')(s^k,\tilde{t})} - \sum _{(j,t') \in {\mathcal {N}}^{+}(s^k,\tilde{t})} y^k_{(s^k,\tilde{t})(j,t')} \right] \nonumber \\&\quad \qquad \qquad \qquad \forall \ k \in {\mathcal {K}}_{in}, \ \forall \ s^k \in {\mathcal {S}}^k_{in}, \ \forall \ t \ge 1, \end{aligned}$$

(33)

$$\begin{aligned} c_{s_l^k} - \sigma ^t_{s_l^k} \le c_{s_l^k} \, \left[ 1-\alpha (s^k_{l+1},t) \right]&\qquad \forall \, k \in {\mathcal {K}}_{in}, \ \forall \, t = 0, \dots , T,\nonumber \\&\qquad \forall \, s^k_l \in \mathcal {S}^k_{in}, \ \forall \, l = 1, \dots , |\mathcal {S}^k_{in}|-1, \end{aligned}$$

(34)

$$\begin{aligned} \sum _{(j,t') \in {\mathcal {N}}^{-}(s^k,t)} y^k_{(j,t')(s^k,t)} \le c_{s^k} \, \alpha (s^k,t)&\qquad \forall \ k \in {\mathcal {K}}_{in}, \ \forall \ s^k \in {\mathcal {S}}^k_{in},\nonumber \\&\qquad \forall \ t \ge 0, \end{aligned}$$

(35)

$$\begin{aligned}&\rho ^t_{s^k} = \sum _{\tilde{t}=0}^{t} \left[ \sum _{(j,t') \in {\mathcal {N}}^{+}(s^k,\tilde{t})} y^k_{(s^k,\tilde{t})(j,t')}- \sum _{(j,t') \in \mathcal {N}^{-}(s^k,\tilde{t})} y^k_{(j,t')(s^k,\tilde{t})}\right] \nonumber \\&\quad \qquad \qquad \forall \ k \in {\mathcal {K}}_{out}, \ \forall \ s^k \in \mathcal {S}^k_{out}, \ \forall \ t \ge 1, \end{aligned}$$

(36)

$$\begin{aligned} c_{s^k_{l}} - \rho ^t_{s^k_l} \le c_{s^k_l} \, (1-\beta (s^k_{l+1},t))&\quad \quad \forall \, k \in {\mathcal {K}}_{out}, \ \forall \, t \ge 0,\nonumber \\&\quad \quad \forall \, s^k_l \in {\mathcal {S}}^k_{out}, \ \forall \, l = 1, \dots , |\mathcal {S}^k_{out}|-1, \end{aligned}$$

(37)

$$\begin{aligned} \sum _{(j,t') \in {\mathcal {N}}^{+}(s^k,t)} y^k_{(s^k,t)(j,t')} \le c_{s} \, \beta (s^k,t)&\quad \forall \, k \in {\mathcal {K}}_{out}, \ \forall \, s^k \in {\mathcal {S}}^k_{out},\nonumber \\&\quad \forall \, t \ge 0. \end{aligned}$$

(38)

Recalling that vehicles F1 and F2 may only move in their respective subgraphs, constraints (18) and (19) ensure they start their route from their parking areas (i.e., $\omega ^{1}$ or $\omega ^{2}$, respectively) at the beginning of the time horizon (i.e., at $t=0$), and return there at the end of the time horizon (i.e., at $t=T$). Constraints (20)–(21) impose that at most one vehicle, either of F1 or F2, can be present in any arc of their respective subgraph. The only exceptions are for the holding arcs representing dwell time at their respective parking areas and, for F2 only, also for the charging station. Note that, compared to the model described in Lanza et al. (2022d), these constraints have been here modified to allow more than one vehicle to idle at the charging station.

Constraints (22) are the flow conservation constraints for the incoming product types $k \in {\mathcal {K}}_{in}$. New releases during the time horizon are represented by values $d_{in}^k(r,t) > 0$ for some r and time instant t. For $t=0$, it is considered the chance of already having some items idling on some r or some b, as a result of operations previously performed. Moreover, notice that the flow of items of a product type $k \in {\mathcal {K}}_{in}$ always terminates in one of its preassigned storage locations. Relations (23) are the flow conservation constraints for $k \in {\mathcal {K}}_{out}$. Here, it is considered the chance of having some items of product type $k \in {\mathcal {K}}_{out}$ idling on some b at time $t=0$, as a result of operations previously performed. Moreover, the flow of items of a product type $k \in {\mathcal {K}}_{out}$ can never be stored in any storage location once retrieved and thus terminates in the collection area. Constraints (24) ensure that all the items of product type $k \in {\mathcal {K}}_{out}$ requested at time t are transported to the collection area before t. Constraints (25) instead define the composition of the collection area at the beginning of the time horizon.

Linking capacity constraints for vehicles of F1 and F2, respectively, are defined by relations (26)–(27). They state that freight flows can only be transported using vehicles which have been selected to move within the warehouse and that the total commodity flow on any moving arc cannot exceed the capacity of the vehicle travelling along it.

Relations (28)–(32) define the capacity constraints for each location of the warehouse. Specifically, (28), (29), (30) relate to input points, to collectors and the collection area, respectively. Moreover, constraints (31) impose that the number of items of product type $k \in {\mathcal {K}}_{in}$ that can be stored into a storage location preassigned to k does not exceed its capacity, while constraints (32) guarantee that the number of items of product type $k \in {\mathcal {K}}_{out}$ that can be retrieved from a storage location occupied by product type k does not exceed the number of items stocked in it.

Constraints (33)–(35) and (36)–(38) regulate the specific storage/picking policy applied in the warehouse prescribing that a new storage location, say $s^k_{l+1}$, may be utilized for storing/picking only if the previous one, say $s^k_{l}$, in the considered order of precedence is already full/empty. For the former, equations (33) define the total number of items of product type $k \in {\mathcal {K}}_{in}$ stocked in the storage location $s^k \in {\mathcal {S}}^k_{in}$ until time t (note that, at $t=0$ and for the first storage location in the given order of precedence, this is an input data). If storage location $s^k_{l}$ has not already reached its saturation at time t, constraints (34) do not allow the next assigned storage location in the related order of precedence, i.e., $s^k_{l+1}$, to be used to stock items of product type $k \in {\mathcal {K}}_{in}$: this is mathematically guaranteed by forcing $\alpha (s^k_{l+1},t) = 0$ in this scenario thanks to constraints (34). As opposed, when storage location $s^k_{l}$ has reached its saturation, i.e., $c_{s_l^k} = \sigma ^t_{s_l^k}$, storage location $s^k_{l+1}$ becomes eligible to stock items of product type k, being $\alpha (s^k_{l+1},t)$ allowed by the combination of constraints (34) and (35) to assume value 1, that is $\alpha (s^k_{l+1},t) = 1$. For the latter, the logic is similar. Equations (36) define the total number of items of product type $k \in {\mathcal {K}}_{out}$ retrieved from the storage location $s^k \in {\mathcal {S}}^k_{out}$ until time t (also in this case, at $t=0$ and for the first storage location in the given order of precedence, this is an input data). Constraints (37) impose that the next storage location in the related order of precedence, $s^k_{l+1}$, cannot be used to retrieve items of product type k, unless the previous storage location, $s^k_{l}$, has been completely emptied. In the latter case, $\beta (s^k_{l+1},t) = 1$ is allowed by the combination of constraints (37) and (38); otherwise $\beta (s^k_{l+1}, t) =0$ and retrieval has still to be performed from $s^k_{l}$.