This paper considers the problem of scheduling automated guided vehicles (AGVs) with battery constraints. Each transport request involves a soft time window, and the AGV fleet used to service those requests is heterogeneous with a diverse set of capabilities and travel costs. In contrast to the existing literature, each transport request may require different AGV material handling capabilities (such as lift loads, tow loads, or handle loads with a mounted robot arm), and the AGV batteries can be recharged partially under consideration of a critical battery threshold. The problem is to assign the transport and charging requests to AGVs, sequence the requests, and determine their starting times and the recharging durations of the AGVs with the objective of minimizing a weighted sum of the tardiness costs of transport requests and travel costs of AGVs. A mixed-integer linear programming model is formulated.
We also propose a new matheuristic that makes use of an adaptive large neighborhood search algorithm and a linear program to solve industry-size instances. We illustrate the efficacy of our approach with an industry case study using real-world data.
1 Introduction
A long-term goal of smart industry is to have a lights-out factory, which requires the control of operations at the shop floor without human interference with the use of automated guided vehicles (AGVs). The AGVs are driverless vehicles that transport goods and materials throughout various areas, e.g., shipping and receiving areas, storage, and workstations ( Confessore, Fabiano, & Liotta, 2013; Vivaldini, Rocha, Martarelli, Becker, & Moreira, 2016 ).
The applications of AGVs have grown enormously because of their benefits, such as flexibility in processes, space utilization, product safety, and computer integration and control. More recently, the AGV technology has enhanced with new advances in their flexibility, where each vehicle will have specific capabilities (e.g., one can lift light or heavy loads, while another can tow loads) to perform heterogeneous tasks ( De Ryck, Versteyhe, & Debrouwere, 2020; Riazi, Diding, Falkman, Bengtsson, & Lennartson, 2019 ). This type of tasks is prevalent in the high-tech manufacturing sector, and it necessitates the handling of materials with specialized equipment, thus requiring a fleet of AGVs with varying capabilities.
An example of a factory requiring different material handling tasks is Brainport Industries Campus (BIC), a new high-tech campus constructed in Eindhoven, The Netherlands, which serves as a home front for far-reaching partnerships between suppliers, specialist companies, and innovative educational and knowledge institutions. The campus is a joint initiative of high-tech suppliers, also called tenants, to co-exist and collaborate on multiple fronts ( Brainport Industries Campus, 2020 ). Since multiple tenants collectively make use of a pool of AGVs, practitioners are looking for a better understanding on how to manage large AGV fleets with heterogeneous capabilities that can cater to a more diverse set of requirements. Consequently, there is an increasing need for planning algorithms to manage such fleets. Industry feedback indicates that most companies rely on the prepackaged fleet management software that is provided by AGV manufacturers. However, such a software is only applicable for a homogeneous AGV fleet and does not make use of look-ahead information such as due dates and time windows. In addition, when considering long travel distances in BIC, the performance of the AGVs can be influenced by battery restrictions. Therefore, efficient management and control of such an AGV system become key factors in reducing material handling costs, work-in-process inventories, and overall operational costs.
In this paper, we address the problem of scheduling AGVs with battery constraints and heterogeneous capabilities. Our problem holds some similarities to the electric vehicle routing problem with time windows (EVRPTW) addressed by Keskin & Çatay (2016) and Zhao & Lu (2019) . The EVRPTW considers a fleet of electric vehicles (EVs) having limited driving range due to their battery capacities, which may need to visit charging stations while servicing customers in their routes. Our problem, however, distinguishes from the EVRPTW literature by the following features to comply with industrial needs. First, the fleet of AGVs considered in this study is heterogeneous in terms of not only travel speeds and costs, charging rates, and discharging rates but also, more importantly, capabilities to service different types of requests from various tenants. Second, partial charging for the AGVs is allowed under consideration of a critical battery threshold. Third, we consider both travel costs of AGVs and tardiness costs of transport requests, where different types of AGVs and requests have different unit travel costs and penalty charges, respectively. Motivated from our industry collaboration at BIC, the combination of these features constitutes the main novelty of the problem. We aim to make decisions on which AGVs the transport requests are assigned, the sequence in which the requests are processed, the times when AGVs should start servicing the requests, and the recharging duration of AGVs at charging stations so that a weighted average of the total travel and total tardiness costs is minimized.
The main contributions of this paper are as follows.
(i) We formulate the problem as a mixed-integer linear programming (MILP) model.
(ii) We propose a new matheuristic approach that surpasses the MILP and the dispatching policy currently used in practice (see Appendix D for details) given the same computational time limit. In our matheuristic, which combines adaptive large neighborhood search (ALNS) ( Ropke & Pisinger, 2006 ) and linear programming (LP), we propose a new ALNS-based algorithm by introducing two pairs of destroy and repair operator sets to perform exploration and exploitation in distinct phases of the search. Here, we introduce problem-specific heuristics concentrating on the removal and insertion of charging tasks in which a possibility of partial charging for AGVs is allowable.
We choose the ALNS framework to solve our problem since it provides flexibility to hybridize with other approaches to make more efficient implementations for optimization problems. Also, it has recently been demonstrated that ALNS-based matheuristics yield promising results on electric vehicle routing problems ( Keskin, Laporte, & Çatay, 2019; Keskin & Çatay, 2018; Koc, Jabali, & Laporte, 2018 ) owing to the ability to find the best tradeoff between the efficacy of exact approaches and the efficiency of metaheuristics ( Guimares, Klabjan, & Almada-Lobo, 2013 ).
(iii) We provide managerial insights by performing a sensitivity analysis on a real-world case study, which allows us to determine potential room for improvement in practice. In this analysis, we vary the scheduling horizon, tightness of time window, the number of requests, AGV fleet size, costs related to tardiness and AGV travel, battery thresholds, and computational time limit. Our computational results reveal that the proposed matheuristic allows finding solutions of high quality even for industry-size large instances and yields an average of 33% improvement with respect to the dispatching policy in practice.
(iv) We provide data for the real-world case study used for the experiments in this paper. The problem instances in the case study are generated from a fieldlab that has been built in the BIC. Additionally, the research is conducted in close collaboration with our industrial partner, IJssel Technology, a company located in the Netherlands, which is the central fleet owner at the BIC. The outcomes of this research can be broadly applied to settings when there is a central operator of a heterogeneous AGV fleet servicing different types of transport requests.
The remainder of this paper is organized as follows. Relevant literature is discussed in Section 2 . The problem is formally described in Section 3 , and a mathematical model is formulated in Section 4 . A matheuristic is proposed in Section 5 to solve industry-sized instances. Next, we describe the industry case study, conduct the parameter tuning, evaluate impacts of operators, and present computational results in Section 6 . Finally, conclusions and future research directions are discussed in Section 7 .
2 Literature review
The design of an AGV control and management system requires the best strategies to solve problems associated with common functions such as dispatching, routing, and scheduling ( Qiu, Hsu, Huang, & Wang, 2002 ). Langevin, Lauzon, & Riopel (1996) define dispatching as the process of selecting and assigning tasks to vehicles, routing as the selection of the specific paths that each vehicle traverses to accomplish its transportation tasks, and scheduling as the determination of the arrival and departure times of vehicles at each station. The decisions related to these functions can be made either simultaneously or sequentially. The reviews of Le-Anh & De Koster (2006) and De Ryck et al. (2020) present guidelines for the control of AGVs and include often-neglected areas, such as idle-vehicle positioning and battery management.
Most dispatching rules in the literature are single-attribute ones, which dispatch vehicles based on only one parameter ( Le-Anh, van der Meer et al., 2004; Li & Kuhl, 2017 ). Commonly used dispatching rules are first-come-first-serve (FCFS) and shortest-travel-distance-first (STTD). While single-attribute dispatching rules are common in practice, multi-attribute dispatching rules, which aggregate different attributes to select a task to be serviced by a vehicle, prove to be better in general ( Le-Anh & De Koster, 2005 ). The methodology proposed in Ho, Liu, & Yih (2012) involves concurrently solving the pickup-dispatching problem along with the load-selection problem with a scoring function that makes use of the slack time of parts, waiting time of parts, and the distance to the load. Jeong & Randhawa (2001) define a scoring function governed by a trained neural net to calculate priorities of incoming jobs. This method is particularly useful if there is a lot of data available on jobs and no clearly defined score function to assess priorities.
Dispatching rules have proven to be good and easy to implement but have been criticized for having a myopic view.
More advanced approaches have, therefore, been studied to improve system performance. Liu, Tan, Kurniawan, Zhang, & Sun (2018) develop a mixed-integer programming model that solves the scheduling problem of mobile robots, tasked with transporting materials among different places, in a periodic or event-driven or hybrid way (having a static schedule being modified with new events). Sabuncuoglu & Kizilisik (2003) study reactive scheduling problems in a flexible manufacturing environment, where the arrival of jobs is not known in advance. They develop a simulation-based scheduling approach and compare offline and online scheduling schemes. Ebben, van der Heijden, & van Harten (2005) present a serial scheduling method for a dynamic resourceconstrained vehicle scheduling problem and analyze its dynamic behavior using discrete event simulation. Nevertheless, these studies do not consider time windows of transport requests and a heterogeneous fleet of AGVs with battery management.
Another stream of relevant work concerns the electric vehicle routing problem (EVRP). It is a variant of the vehicle routing problem (VRP) in which a fleet of EVs is used, instead of internal combustion engine vehicles (ICEVs). Conrad & Figliozzi (2011) study the EVRP, where EVs may recharge fully or 80% of the battery capacity at certain customer locations during their services. The objectives are to minimize the number of vehicles and total costs related to travel distance, service time, and charging. They propose an iterative construction and improvement heuristic to solve the problem. Erdoan & Miller-Hooks (2012) then study the green VRP (GVRP) in which alternative fuel vehicles such as solar, electric, or biodiesel vehicles are fully refueled at stations en route within a constant amount of time. They propose two heuristics to solve the GVRP with the objective of minimizing the total travel distance.
Schneider, Stenger, & Goeke (2014) introduce the EVRPTW where the recharging may take place at any battery level with a linear charging function, and EVs always leave stations with a full battery.
In their work, a hybrid heuristic of variable neighborhood search and tabu search (TS) is proposed for solving the EVRPTW given the objective being to minimize the total travel distance by using the minimum number of vehicles.
Recently, many authors have relaxed the full charge restriction and allowed partial recharging ( Bruglieri, Mancini, Pezzella, Pisacane, & Suraci, 2017; Bruglieri, Pezzella, Pisacane, & Suraci, 2015; Ding, Batta, & Kwon, 2015; Keskin & Çatay, 2016; Roberti & Wen, 2016; Schiffer & Walther, 2018b; Wen, Linde, Ropke, Mirchandani, & Larsen, 2016 ) with different charging technologies at stations ( Felipe, Ortuño, Righini, & Tirado, 2014; Keskin & Çatay, 2018 ). Several other extensions have also been studied. Some models take into account decisions on both vehicle routing and locating charging stations ( Schiffer & Walther, 2017; 2018b ). Other studies have considered the possibility of a heterogeneous fleet.
Some of them analyze a fleet of EVs having different characteristics such as battery capacities and operating costs ( Hiermann, Puchinger, Ropke, & Hartl, 2016; Zhao & Lu, 2019 ), while others deal with a mixed fleet of EVs and ICEVs ( Goeke & Schneider, 2015; Kopfer & Vornhusen, 2019; Sassi, Cherif-Khettaf, & Oulamara, 2015a; 2015b ). Additionally, several works look at penalties for late arrivals ( Grandinetti, Guerriero, Pezzella, & Pisacane, 2016; Keskin et al., 2019 ). Nevertheless, to the best of our knowledge, the problem in this paper, which concurrently considers transport requests with soft delivery times, heterogeneous AGVs with different material handling capabilities, and partial charging with a critical battery threshold, has not been addressed in the literature. We describe the considered problem in detail in the next section.
R set of transport requests
B set of charging tasks (requests)
T set of all transport requests R and charging tasks B ,where T = R ∪ B
ARr set of capability requirements of request r
SRr pickup or source node for request r
DRr drop-off or destination node for request r
cRr penalty cost incurred per time unit of tardiness for request r
eRr earliest pickup time of a request r
lRr latest delivery time of request r
gRr type of task (request) r
V set of automated guided vehicles (AGVs)
Avk set of capabilities of AGV k
cvk travel cost per time unit of AGV k
charging rate of AGV k
cvk discharging rate of AGV k
bvk critical battery threshold of AGV k
bvk non-critical battery threshold of AGV k
bvk charge level of AGV k
Bvk initial charge level of AGV k
bl minimum battery charge level (%)
bu maximum battery charge level (%)
N set of all nodes in a layout
E set of charging nodes in a layout
hni service time at a node i
dij distance between node i and node j
tijk time to travel between node i and node j by vehicle k
Qk sequence of tasks of AGV k
Q set of task sequences of k AGVs, where Q ={ Q 1 , Q 2 , . . . , Q k }
Sk schedule of tasks of AGV k
S set of task schedules of k AGVs, where S = { S 1 , S 2 , . . . , S k }
Kr set of AGVs having the required capability to process request r
a weight factor
q number of requests to be removed
p degree of randomness in choosing related tasks
k level in the class of regret- κ heuristics
λ reaction factor to changes in the effectiveness of operators
w initial temperature control parameter
ε cooling rate in the simulated annealing
φ1 weight for relatedness in terms of distance
φ2 weight for relatedness in terms of time-window
φ3 weight for relatedness in terms of AGV capability requirement
ψ1 score of a heuristic rewarded when finding a new best solution
ψ2 score of a heuristic rewarded when finding a better solution
ψ3 score of a heuristic rewarded when finding a worse solution but accepted
nc number of iterations station methods run when in- voked
Tmax execution time limit of the matheuristic
γ number of random tasks to remove
3 Problem description
The problem in this paper concerns a set of transport requests R that is serviced by a heterogeneous AGV fleet V with battery constraints. The source (pickup) and destination (delivery) nodes of reR quest r ∈ R are denoted with s R r and d r , respectively. Each request r ∈ R requires a set of capabilities A R r from an AGV for transport, R ], where e R is the earliest pickup and it has a time window [ e R , l r r r time and l r R is the latest delivery time of the request. An AGV is not allowed to reach the source of a request before its earliest pickup time. On the other hand, if AGV k arrives at the destination of request r after its latest delivery time, the tardiness occurs and is measured by max { 0 , (arrival time of AGV k at d r R ) − l r R } , with a penalty cost c r R incurred per time unit.
Each request needs to be performed by exactly one capable AGV in the heterogeneous fleet. The AGV k ∈ V has a set of capabiliV ties A V k and is capable of servicing request r ∈ R if A R r ⊆ A k holds, i.e., the capability requirements of a request are included in the capabilities of an AGV. For example, with a request classified as a heavy load that needs to be lifted on a pallet, only those AGVs that have the capability to lift heavy loads can carry out that request.
Each AGV also has a travel cost per time unit c k V . While an AGV is traveling, the battery charge level decreases proportionally with the traversed distance at a discharging rate of d ˆ k V , and the AGV may need to visit a charging station before continuing operations.
Here, the battery is recharged at a charging rate of c ˆ k V . The charge level of an AGV (measured as a percentage) must take values between the minimum and maximum battery charge level denoted by b l percent and b u percent, respectively. The AGV must recharge its battery if the charge level drops below a critical threshold b k , and the recharging duration should be long enough to allow the charge level to reach at least this threshold. In other words, before starting a new request, the charge level of the AGV must V be at or above b k percent. The critical threshold is a prespecified parameter for a given network of nodes to assure that an AGV can reach any request destination and a charging station right afterwards if needed. It is notable that an AGV is not permitted to visit a charging station while performing a request (i.e., carrying a load), and its route starts and ends at one of the charging stations. An AGV also cannot service more than one request at the same time, i.e., on its tour from a source node to a destination node (singleload AGVs). In addition, the AGVs can detect objects on the shop floor and move around them while traveling. Consequently, collisions between AGVs can be avoided by hardware, thus not being considered in this paper.
The problem can be defined on a directed graph G = ( N, A ) , where N = { 1 , 2 , . . . , n } is the set of nodes and A = { ( i, j ) | i, j ∈ N, i = j } is the set of arcs. A node in the graph can represent a pickup or delivery location, or a charging station. Each pickup or delivery node i ∈ N has a service time h N that represents the time i needed for loading or unloading activities, while each charging node has zero service time. Each arc ( i, j ) is associated with a distance d i j , and it is traveled by AGV k ∈ V in t i jk time units. Note that both distance and travel time are symmetrical, i.e., d i j = d ji and t i jk = t jik . In addition, let E ⊂ N denote the set of charging stations and B = { 1 , 2 , . . . , n B } denote the set of charging requests, where n B is a safe upper bound, computed as n B = | V | · | E | , which allows AGVs to charge as many times as needed. Therefore, not all the requests from B need to be used in a schedule. For each charging request r ∈ B , the source and destination nodes are the same charging station, and the earliest pickup and latest delivery time R are set to zero and infinity, respectively ( e R r = 0 , l r = ∞ , ∀ r ∈ B ).
Lastly, T denotes the set of all transport and charging requests such that T = R ∪ B . Note that R , V, and N in the superscript are referred to request, AGV, and node, respectively.
The goal of the problem is to determine (a) an assignment of AGVs to requests including charging requests, (b) a sequence of requests on AGVs, and © starting times of requests and recharging duration such that problem constraints are satisfied, and a weighted sum of the total tardiness cost of requests and the total travel cost of AGVs is minimized.
4 Mathematical formulation
In this section, we formulate a mathematical model based on the description in Section 3 . The model formulation is presented below.
4.1. Decision variables
In addition to the notation presented in Section 3 , we introduce the following decision variables for the mathematical formulation:
x r r k binary variable, equal to 1 if request r is performed prior to r by AGV k , 0 otherwise
y S rk arrival time of AGV k at the source node of request r ( S refers to source)
y Drk arrival time of AGV k at the destination node of request r ( D refers to destination)
T rk tardiness of request r performed by AGV k
z rk percent amount of battery discharge of AGV k when it reaches the source node of request r
δ rk travel time of AGV k when it performs request r
4.2. Mixed-integer linear programming model
The mathematical model of the described problem can be formulated as follows:
The objective function (1) minimizes the weighted sum of tardiness costs of requests and travel costs of AGVs, where α is a weight coefficient to prioritize one objective over the other. Constraints (2) ensure that each request is attended by at most one AGV, while Constraints (3) define whether or not a charging request is performed. Constraints (4) ensure that the capabilities of the requests and AGVs match. Self-visits are avoided by Constraints (5) .
Constraints (6) conserve the incoming and outgoing arcs for each request done by an AGV. Constraints (7) ensure that an AGV arrives at the source node after the earliest pickup time. Constraints (8) and (9) compute the arrival times at destination and source nodes, respectively. Also, Constraints (9) prevent subtours with M defined as a large positive constant. Constraints (10) and Constraints (18) ( T rk ≥ 0 ) define the tardiness of each request. Constraints (11) determine the travel time of an AGV to perform a request. Constraints (12) and (13) calculate the amount of battery discharge due to the travels between the source and destination of a request and between the destination and source of two requests, respectively. The charged amount due to a charging request is given by Constraints (14) . Constraints ( 15–17 ) set the lower and upper limits ( b l = 0 , b u = 100 ) for an amount of battery discharge, while Constraints ( 18 –19 ) guarantee valid domains for the remaining decision variables.
5. Proposed matheuristic
In this section, we propose a matheuristic based on the integration of ALNS and LP to solve the problem in Section 3 . In our matheuristic, the ALNS is responsible for solving the assignments of AGVs and sequencing of requests, while the LP model is used to determine the schedule of requests, i.e., starting times and recharging duration. The ALNS involves multiple destroy and repair operators, which compete to modify a current solution. In each iteration, a destroy and a repair operator are chosen according to their past performance, and the resulting new solution is accepted if it satisfies some criteria defined by the ALNS ( Demir, Bekta, & Laporte, 2012 ). Here, we integrate the LP into the ALNS framework to utilize the advantages of both methods to find high-quality solutions for large problem instances efficiently.
Let Q 0 and S 0 denote the initial sequence and initial schedule, respectively. Similarly, Q current and S current denote the current sequence and schedule, while Q best and S best denote the best sequence and schedule, respectively. In the proposed matheuristic, there are four sets of operators: request destroy, request repair, (charge) station destroy and station repair, which are denoted by + − + − + − c , c , s , and s , respectively. Let d c ∈ c and r c ∈ c denote a destroy and a repair operator for requests, respectively, while + d s ∈ − s and r s ∈ s denote a destroy and a repair operator for charging tasks, respectively. Also, the weights of operators in − c , − + + c , s , and s , which reflect probabilities to choose an operator, + − + are stored in vectors P − c , P c , P s , and P s , respectively. The pseudocode of the proposed matheuristic is presented in Algorithm 1 .
The matheuristic starts with the initialization of a feasible se+ − + quence Q 0 and weights of operators in P − c , P c , P s , and P s (line 1).
Initially, all operators have the same weight. The initial schedule S 0 is then computed by an LP model with the input of Q 0 (line 2). The matheuristic considers Q 0 ( S 0 ) the current and best so far (lines 3– 4) and sets the simulated annealing (SA) temperature τ to its initial value τ 0 (line 5). Then, a loop is executed until a stopping criterion is met (lines 6–31). At the start of the loop, the matheuristic selects a pair of request destroy and repair operators d c and r c to modify the current sequence Q current to generate a new sequence Q new , which is next used to compute S new (lines 7–9). If a new best solution is found, the matheuristic focuses on exploiting its neighborhood by applying station destroy and repair operators d s and r s until it reaches n c iterations and does not obtain another new best solution (lines 10–25). Otherwise, a new solution is accepted if it is better than the current solution, or worsen than but satisfies the SA acceptance criterion, where r U ( 0 , 1 ) is a random number generated according to the uniform distribution U ( 0 , 1 ) (lines 26– 27). The weights of selected request and station operators in P − c , − + P + c and P s , P s are updated when being used, respectively (lines 23 and 29). Also, the SA temperature is the last to be updated in each iteration, where ε is the cooling rate (line 30). At the end of the matheuristic, the best schedule S best is returned (line 32). We provide details on the algorithm in separate sections, as stated in the pseudocode
5.1. Solution representation
The solution of our problem is a schedule in which tasks are carried out at specific time moments. In this section, we describe a solution representation by presenting the types, sequence, and schedule of tasks subsequently as follows.
5.1.1. Types of tasks
A schedule S consists of three types of tasks such as transport request (loaded travel) T R , unloaded travel UT , and charging CH, which are depicted in Fig. 1 . Let g R r denote the type of task r.
First, the task in Fig. 1 (a) represents request r of type T R ( g R r = T R ), R which has the time window [ e R r , l r ], source and destination nodes R R ( s R r , d r ), and set of capability requirements A r . Second, Fig. 1 (b) presents a task of type UT ( g R = UT ), which indicates the travel r of an AGV without any load to move from the destination of its previous task ( s R r ) to the source of its next task ( d r ). This travel occurs in between two tasks of type T R or types T R and CH. Finally, Fig. 1 © presents a task of type CH ( g R r = CH), which has the source and destination nodes at the same charging station, i.e., R R R s R r = d r , where s r , d r ∈ E.
5.1.2. Sequence of tasks
Task sequences on AGVs are constructed from tasks of two types T R and CH. We denote by Q the set of task sequences of all AGVs, i.e., Q = { Q k | ∀ k ∈ V } , where Q k is the sequence of tasks on AGV k . Additionally, let Q r̄ k denote the r̄ th task in sequence Q k , and the attributes of the task are represented by ( e R , l r̄ R , s R , d r̄ R , A R , g R ) . For example, Fig. 2 illustrates the set Q r̄ r̄ r̄ r̄ consisting of two AGVs, with their set of capabilities, and their respective task sequences Q 1 and Q 2 in which Q 3 1 and Q 2 2 are the third and second task of AGV 1 and 2, respectively. Here, the attribute values of Q 3 1 are ( −, −, 19 , 19 , −, CH ) , while those of Q 2 2 are ( 32 , 55 , 1 , 12 , { A, C } , T R ) . Note that the symbol “-” indicates the absence of the corresponding attribute of a task.
5.1.3. Schedule of tasks
To compute the schedule S of all AGVs, each sequence Q k is post-processed by adding tasks of type UT to result in sequence Q ˆ k such that it is composed of transport requests, charging tasks, and unloaded travel tasks. Let S k denote the schedule of tasks on AGV k , then S = { S k | ∀ k ∈ V } . We also denote by S r̄ k the r̄ th task in schedule S k . Two additional attributes, namely y S r̄ and y D , are r̄ included to indicate the arrival times of an AGV at the source and destination nodes of task r̄ , respectively. Note that these time moments are determined while computing the schedule S (see Section 5.6 ). In addition, for a task of type CH, the difference between y D and y S r̄ (i.e., y D − y S r̄ ) implies the recharging durar̄ r̄ tion. Fig. 3 illustrates the set S = { S 1 , S 2 } , where S 1 and S 2 are the schedule of tasks on AGV 1 and 2, respectively. In schedule S 1 , S 2 1 indicates the second task ( r̄ = 2 ) with its attributes ( e R , l r̄ R , s R , d r̄ R , A R , g R , y S r̄ , y D ) = ( 8 , 25 , 5 , 11 , { A, B } , T R, 8 , 19 ) .
r̄ r̄ r̄ r̄ r̄ After obtaining the schedule S, its objective value is calculated according to Eq. (20) , which is equal to the weighted sum of tardiness costs of requests and travel costs of AGVs. Here, ( y D − l r̄ R ) + r̄ gives the tardiness of each transport request r̄ .
formula 20
5.2. Initialization
The quality of the initial sequence may have a crucial impact on the final outcome of the matheuristic. To generate a good initial sequence, we propose a constructive heuristic, which employs the Earliest Due Date (EDD) rule and the critical charge insertion method described in Section 5.5.2 . The pseudocode of this constructive heuristic is given in Algorithm 2 . It begins with sorting all transport requests r ∈ R in non-decreasing order of their latest delivery times to create the sorted set R ˆ (line 1). Then, we create the set of capable AGVs, denoted by V C , for servicing each sorted request (lines 2–3). If there exists only one candidate in V C , it is selected (lines 4–5). Otherwise, the AGV that has the minimum additional cost when inserting request r to its sequence is chosen (lines 6–10). The cost function f ( r, k ) of inserting request r to AGV k is determined in line 8, where the tentative arrival time y ˆ D r of the AGV at the request destination is calculated given that the AGV should reach the source node as early as possible but not before the earliest pickup time. After servicing the request, if the charge level b V v of the chosen AGV v drops below the critical battery threshold b v , a charging task is assigned to that AGV.
The (re)charging duration should allow the charge level to reach at least the critical threshold, as mentioned in Section 3 . In general, the proposed constructive heuristic can be classified as a greedy sequencing algorithm.
5.3. Adaptive weights for operators
5.4. Destroy operators
In this section, we present nine destroy operators used in our matheuristic. The first five are classified as request destroy that removes transport requests, whereas the last two are classified as station destroy that removes charging tasks. These operators are either adapted or inspired by Ropke & Pisinger (2006) , Keskin & Çatay (2016) , and Shaw (1998) . In general, the destroy phase mainly consists of removing a number of requests from the current sequence, adding them into a removal list, and returning a partially destroyed sequence.
5.4.1. Request destroy operators
We now describe the operators removing requests as follows.
Shaw destroy (SHD) : This operator is designed to remove requests which are related to each other in respect of several aspects and thus are easy to reshuffle them around. In our algorithm, the relatedness between two requests r and r is defined by S ( r, r ) and is computed by Eq. (23) :
formula23
where φ 1 , φ 2 , and φ 3 are weights (Shaw parameters) corresponding to the distance, time window, and capability requirement terms, respectively, and K r represents the set of AGVs that are capable of servicing request r. A lower S ( r, r ) indicates a higher degree of relatedness between two requests r and r . It is assumed R that d xy , e R x , and l x are normalized by scaling those such that they only take the values between [0,1], where x and y respectively represent the first and second nodes in the distance term, while x in the time-window term indicates either request r or r . The SHD first chooses randomly a transport request r to be removed from a sequence (set) Q. Then, it sorts an array of non-removed requests L in increasing order of their relatedness values with respect to request r. Next, it selects the request placed in the position defined by [ y p | L | ] , where y ∈ [0 , 1 ) is a random number, and p is a deterministic power introducing some variation in the selection of related requests, making related requests more likely to be selected.
[ ·] is rounded to the nearest integer. The procedure iterates in a similar fashion until q requests are removed ( Keskin & Çatay, 2016; Ropke & Pisinger, 2006 ). The value of q will be determined in the parameter tuning section ( Section 6.2 ).
We also use three special cases of the SHD operator referred to as distance- based Shaw destroy (SHD-D), time-window- based Shaw destroy (SHD-T), and capability- based Shaw destroy (SHD-C), where φ 1 = 1 , φ 2 = 1 , and φ 3 = 1 , respectively, and the other parameters in Eq. (23) are set to 0.
Worst request destroy (WRD) : This operator attempts to remove q requests with a high cost from a sequence Q. We define the cost R + of a request r in a schedule S as f ( r, S ) = ( y D r − l r ) . The operator reuses some of the ideas from the SHD. At each iteration, it sorts an array of non-removed requests L in decreasing order of their costs and removes the request listed in position [ y p | L | ] , where y and p are defined as in the SHD. This removal process is repeated for q iterations.
Random request destroy (RRD) : The operator simply removes γ requests randomly chosen from the sequence Q to diversify the search mechanism.
5.4.2. Station destroy operators
Charging tasks at stations are crucial components of our problem. Hence, removing or repositioning them in an AGV sequence may improve the solution. We here use the following station destroy operators.
Worst charge destroy (WCD) : This operator chooses an AGV k randomly and calculates the amount of charge that it accumulates at each particular charging station. Then, the WCD selects the charging task that leads to the least increase in charge and removes it from the sequence Q k .
Random charge destroy (RCD) : The operator simply selects an AGV randomly and removes a charging task which is also chosen randomly from the sequence of that AGV.
All charge destroy (ACD) : The operator simply selects an AGV randomly and removes all charging tasks from the sequence of that AGV.
5.5. Repair operators
To repair a partially destroyed sequence, six insertion operators are employed in our algorithm. These operators are divided into two groups, namely request repair and station repair. While the operators in the first group attempt to reinsert transport requests removed by destroy methods, those in the second group insert charging tasks into partial sequences thus affecting the visits to the charging stations of AGVs along their route. The insertion operators return a feasible sequence Q that serves as the input for the LP model, presented in Section 5.6 , to compute a feasible schedule S.
5.5.1. Request repair operators
We adapt the greedy and regret- κ insertions from Ropke & Pisinger (2006) and propose a new method called EDD-based greedy insertion to reinsert the transport requests. These three operators are described as follows.
Greedy insertion (GRI) : This operator repeatedly inserts a request into its best position overall that is called the minimum cost pok sition. Let Q r→ denote the resulting sequence when inserting rei quest r at position i in sequence Q k and f ( r, k, i ) denote the cost of the whole partial sequence Q after the insertion. The pseudocode of the greedy insertion is presented in Algorithm 3 . For each request r in the removal list D , the operator collects a set of capable AGVs V C for servicing that request (line 3) and determines the cost f ( r, k, i ) (line 5), where y ˆ D is the tentative arrival time at the rer̄ quest destination as defined in Section 5.2 . The process is repeated for all requests in the set D , and the request that has the minimum cost f ( r, k, i ) is inserted at its designated position (lines 8–9). The process continues until all destroyed requests are reinserted.
Regret- κ insertion (RKI) : One potential downside of the GRI is that it often postpones the insertion of requests that have high insertion costs. When those requests need to be reinserted, there tion, the operator chooses the request that has the maximum regret value overall (line 12), where each regret value f regret− κ ( r ) is determined in line 10. Then, the chosen request is inserted into its minimum cost position. The process is repeated until all removed requests are reinserted. We use small values of κ , i.e., the regret-2 and regret-3 heuristics are employed in our matheuristic.
EDD-based greedy insertion (EGI) : Although the greedy and regret- κ insertions are often the preferred heuristics in the ALNS literature ( Demir et al., 2012; Grangier, Gendreau, Lehd, & Rousseau, 2016; Gullhav, Cordeau, Hvattum, & Nygreen, 2017; Keskin & Çatay, 2016; Žulj, Kramer, & Schneider, 2018 ), they are computationally expensive since the insertion cost has to be calculated for each available position in a sequence Q. We hence introduce a faster heuristic called the EDD-based greedy insertion. Here, the difference between the EGI and the GRI is that the EGI only considers to reinsert a destroyed request r ∈ D into the sequence of a capable AGV k ∈ V C at the first position i , where l r R < l R k . The idea lies in the fact that requests with earlier delivery times (earlier due dates) should be serviced before the ones with later delivery times.
It is notable that a request insertion operator is feasible with respect to capability requirements but may not meet the critical battery threshold. In this case, the critical charge insertion presented in Section 5.5.2 is applied to make repaired solutions feasible.
5.5.2. Station repair operators
After removing charging tasks, a partially destroyed sequence may become infeasible to the battery charge requirement of AGVs.
To repair the sequence, we propose two operators inserting charging tasks, namely critical charge and non-critical charge insertions.
Since AGVs can charge as many times as needed, any charging task can be inserted throughout the station repair operators, and thus it may not be necessary to reinsert all the charging tasks which are removed by the station destroy heuristics. Note that any AGV should perform a charging task at a charging station closest to the location of its previously serviced request.
Critical charge insertion (CCI) : This operator estimates the battery charge level after servicing each request in a partial sequence
