Task scheduling for transport and pick robots in logistics: a comparative study on constructive heuristics

We study the Transport and Pick Robots Task Scheduling (TPS) problem, in which two teams of specialized robots, transport robots and pick robots, collaborate to execute multi-station order fulfillment tasks in logistic environments. The objective is to plan a collective time-extended task schedule with the minimization of makespan. However, for this recently formulated problem, it is still unclear how to obtain satisfying results efficiently. In this research, we design several constructive heuristics to solve this problem based on the introduced sequence models. Theoretically, we give time complexity analysis or feasibility guarantees of these heuristics; empirically, we evaluate the makespan performance criteria and computation time on designed dataset. Computational results demonstrate that coupled append heuristic works better for the most cases within reasonable computation time. Coupled heuristics work better than decoupled heuristics prominently on instances with relative few pick robot numbers and large work zones. The law of diminishing marginal utility is also observed concerning the overall system performance and different transport-pick robot numbers.


Introduction
In the Heterogeneous Robotic Order Fulfillment System (HROFS), two types of robots with specialized and complementary capabilities exist: 1) transport robots with object transfer and transient storage capability, typically automated guided vehicles (AGVs) or autonomous mobile robots (AMRs), that autonomously receive order fulfillment tasks, travel across the workspace and retrieve items from different storage stations and 2) pick robots with mobile manipulation capability, typically mobile manipulators, mobile dual-arm robots, hybrid leg-wheel robots or even human workers enhanced by mobile platforms (only for tricky products that is still challenging for robots), that fetch items from storage stations and place them into transport robots' containers. As is verified in our previous research [1], by combining two types of robots, practitioners can achieve the same or better system performance with relatively lower cost.
Targeting the task allocation and scheduling aspect of the HROFS, we formulate the Transport and Pick Robots Task Scheduling (TPS) problem, which seeks the collective time-extended task schedule for multiple transport robots and multiple pick robots collaboratively performing order fulfillment tasks. In the TPS problem, as customer orders are characterized by multiple lines of miscellaneous items, both two robot types have to travel to different locations to collect items or perform pick&place operations. Furthermore, as we split object transfer and pick&place operations and assign them to two robot types, to achieve better performance, it is paramount to enable them to act tightly and synergistically. According to the multi-robot task allocation (MRTA) taxonomy iTax [2], the TPS problem falls in the complexschedule [CD] category, with single-task robots [ST], multiple-robot tasks [MR], and the time-extended allocation [TA] problem (CD[ST-MR-TA]). iTax by Korsah, et al. [2] is the first taxonomy for multi-robot task allocation that explicitly takes issues of inter-task constraints and inter-robot utilities into consideration. Nunes, et al. [3] further propose a novel taxonomy for MRTA problems with temporal and ordering constraints (MRTA/TOC), which belongs to the XD[MT-SR-TA] category in iTax. Typical MRTA methods include centralized methods (exact [4], approximate [5], heuristics [6] and metaheuristics [7]), decentralized methods (game-theoretic [8], market and negotiation-based [9]), swarm-based methods [10] and hybrid methods [11,12]. For research on MRTA with complex-schedule constraints, materials are relatively rare. Jones, et al. [13] solve the time-extended multi-robot coordination with intra-path precedence constraints in the disaster response domain. The robot team comprises of multiple unmanned fire trucks and multiple debris cleaning bulldozers. They propose two methods, a hybrid method incorporating tiered auctions and two heuristic techniques, clustering and opportunistic path planning. The other method is genetic algorithm. Strengths and weakness of both methods are analyzed and verified by simulations.
The TPS problem is also an integrated routing and scheduling problem, and several combinatorial problems are related to it. Vehicle routing problem (VRP) [14] and its variants capture the routing components of the TPS problem. Open shop scheduling problem (OSP) [15] and its variants capture the scheduling components of the TPS problem. The open shop scheduling problem also has many integrated routing and scheduling variants, such as routing open shop scheduling [16], open shop scheduling with sequence-dependent setup times [17] and open shop scheduling with transportation/travel times [18]. These problems can be solved by exact methods [19], approximate methods [16,20,21], heuristics [22] and metaheuristics [17,18,23,24].
In our previous work [25] and [1], we use decoupled and coupled methods to plan collective time-extended task schedule respectively. Specially in [1], we propose a rank-minimal heuristic based on tripartite matching. We use genetic algorithm as the solver in each iteration. In practice, decision-making on task level should be made on-the-fly, and obtaining satisfying results in a short time is paramount. For such circumstances heuristic methods are natural choices. In fact, as a typical integrated routing and scheduling problem, there exist many potential heuristics yet to be explored for the TPS problem.
However, there is a lack of a comprehensive heuristic algorithm design principles and it is still unclear how to choose among multiple methods. In view of this, in this research we concentrate on designing constructive heuristics, conduct a comparative study, and figure out the most promising heuristic. This paper is organized as follows. In Section 2, we formulate the TPS problem and introduce the sequence models. In Section 3, we design five constructive heuristics together with complexity analysis. In Section 4, we conduct a comparative study on these heuristics and analyze computational results. Finally, in Section 5, we provide the conclusions and point out possible future work.

Problem formulation
The TPS problem is defined by three sets, a heterogeneous robot set R, a logistic network G and an order fulfillment task set T . We explicitly consider the environment or the logistic network because pick robots are associated with their dedicated work zones.
Logistic Network: The logistic network is represented as a simple, undirected and connected graph G = (V, E). The vertex set V = V P ∪ V D is a union of the pickup station subset V P and the delivery station subset V D , and V P ∪ V D = ∅. A delivery station v ∈ V D is a vertex where transport robots receive order fulfillment tasks, start from and return to. A pickup station v ∈ V P is a vertex where two types of robots meet so that pick robots perform pick&place operations, and transport robots collect items. The edge set E consists of connections between stations and each edge e ∈ E is assigned to a travel cost w(e).
Robots: The heterogeneous robot set R consists of m transport robots and n pick robots, i.e., R = R T i |i = 1, 2, · · · , m ∪ R P j |j = 1, 2, · · · , n . A transport robot R T i starts from a delivery station and undertakes an order fulfillment task from the task set, travels across the whole workspace as subtasks indicates, rendezvous with pick robots to allow them pick&place items, and finally returns to the delivery station after all items are collected. A pick robot R P j is dedicated to a work zone represented as a set of pickup stations Z i , Z i ⊂ V P . It picks items from storage (racks, shelves) at corresponding pickup stations, meets with transport robots, and places items into transport robots' containers. We do not allow any two work zones overlap, i.e., ∀1 ≤ p, q ≤ n, p = q, Z p ∩ Z q = ∅. Meanwhile, all work zones should cover the pickup station subset, i.e., ∪ n j=1 Z j = V P . Tasks: Let T = {T 1 , T 2 , · · · , T m } denote the set of m order fulfillment tasks. Currently for one-shot scheduling we only consider tasks with the same transport robot number m, and we also presume the transport robot R T i has been pre-assigned to the task T i ∈ T . This is because task allocation problems (assigning tasks to transport robots) can be solved independently from task scheduling. Each task T i is a set of several pickup subtasks and one delivery subtask, i.e., T i = {t i00 }∪{t i1k , k = 0, 1, · · · , K i1 }∪ · · · ∪ t ijk , k = 0, 1, · · · , K ij ∪ · · · ∪ {t ink , k = 0, 1, · · · , K in }, where K ij ∈ N + ∪ {0} is the number of pickup subtasks in the zone Z j . A pickup subtask t ijk is defined by its pickup station v ijk and a positive pick&place operate time τ ijk , i.e., Similarly the delivery subtask t i00 is defined by its delivery station v i00 and a positive operate time τ i00 , i.e., t i00 : Pick&place operate time is simply called operate time afterwards. This task definition caters to real-life order-fulfillment tasks in which customer orders may include multiple lines of items, and multiple orders from different customers can also be combined by the order batching process [26].
Outputs: We seek all robots' collective time-extended For subtask t ijk belonging to transport robot R T i 's task, or in pick robot R P j 's work zones, ST * ijk , AT * ijk , PT * ijk , CT * ijk are start time, arrive time, pick time and completion time respectively.

Problem 1 Transport and Pick Robots Task Scheduling (TPS) Problem
Given a TPS problem instance R, G, T as described previously, find a collective timeextended task schedule T S, so that all robots collaborate with each other and accomplish all allocated tasks with the minimization of the makespan C max . If task T i 's completion time is C i , the makespan is C max = max{C i }, i = 1, 2, · · · , m.
A TPS problem instance with 3 transport robots and 4 pick robots is depicted in Fig. 1. In this scenario, 3 transport robots R T 1 , R T 2 , R T 3 and 4 pick robots R P 1 , R P 2 , R P 3 , R P 4 collaborate to execute 3 order fulfillment tasks T 1 , T 2 , T 3 . The objective is to find collective task schedule composed of all 7 robots' schedules, with the makespan minimized.
Assumptions: This research is subject to several assumptions. We concentrate on one-shot deterministic task scheduling with all necessary information available. Travel time estimation is sufficient enough for robots to travel between different subtask stations with collisions or deadlocks avoided. Operate times are adequate for pick robots to perform pick&place operations. Work zones for all pick robots are equal.
According to taxonomy iTax [2], the TPS problem falls in the CD[ST-MR-TA] category. This is because transport robots' and pick robots' task execution orders all are not predetermined; each robot can merely execute one task at the same time; and each subtask requires two types of robots work at the same time. Complex-schedule constraints exist between tasks and the utilities of all robots are interrelated and interdependent. Each robot's schedule is coupled with others and cannot be independently determined. The search space of task execution orders is exponentially large. Two classes of coupledness exist in the TPS problem: coupledness between two types of robots, and coupledness between routing and scheduling components.

Sequence models
To establish the partial relations between two subtasks for all robots unanimously, we can use directed acyclic graph (DAG). As is stated in our previous work [1], the TPS problem and the open shop problem share the same property: two types of orders (transport robot subtask orders and pick robot subtask orders, or job orders and machines orders) should be determined by the task scheduler. For open shop scheduling problems, the permutation list [27] and the rank matrix (can be bijectively mapped to a sequence graph) [20,22,23,[28][29][30][31] are commonly used models. The permutation list is simple and easy to build, however, it suffers from a notorious drawback that is high redundancy: multiple permutation lists possibly mean the same sequence. Oppositely, the rank matrix is non-redundant, yet to guarantee feasibility and nonredundancy, more pre-processing and post-processing algorithms are needed.
For the TPS problem, as the combinatorial space is much more complicated than open shop scheduling problems, non-redundancy property is more preferable. Using the rank matrix models could provide more insightful descriptions on the solution structure. Meanwhile, as they can always provide feasible and non-redundant solutions beforehand, efficient algorithms can be developed. For the TPS problem, block sequence graph and block rank matrix are extended from sequence graph and rank matrix models for open shop problems.
A block sequence graph is feasible if the corresponding directed acyclic graph obtained by the union of all transport robots' subtask orders and pick robots' subtask orders is acyclic. block sequence graph differs from (preemptive) sequence graph [20] in that direct precedence constraints between two split operations need not be considered, as a result, the search space of block sequence graph is larger.
For each block sequence graph, there exists a block rank matrix capturing its structural properties algebraically. In a block rank matrix, an entry, namely rank s ijk means that in the corresponding block sequence graph a path to subtask t ijk with a maximal number of subtasks includes s ijk subtasks. A block rank matrix can always be transformed into a permutation list by topological sorting or linearization. The block sequence graph and its corresponding Fig. 1 A TPS problem instance with 3 transport robots (red, green and blue discs with number 1,2,3), 4 pick robots (white discs with number 1,2,3,4) and 3 tasks (red, green or blue crosses). Transport robots and their corresponding order fulfillment tasks are depicted by the same color (red, green and blue). Pick robots are confined in their corresponding work zones Z 1 , Z 2 , Z 3 , Z 4 , each of which is an aisle between two opposite racks (black bars). Transport robots start from the delivery stations (gray squares in periphery of the workspace), and in contrast, pick robots start from the pickup stations (gray squares near the racks). Each task contains several subtasks possibly in different work zones block rank matrix are shown in Fig. 2 and Eqs. (1), (2).

Constructive heuristics
In this section, we design efficient constructive heuristics for the TPS problem utilizing previous sequence models. Based on whether two types of robots are scheduled at the same time, constructive heuristics can be divided into coupled and decoupled heuristics. Typical discrete operations include matching, append and insertion. Some of these operations have been designed for open shop scheduling problems in [22]. For integrated routing and scheduling problems or variants of open shop scheduling [16,21,32], although progresses have been made to find approximate methods with bounded suboptimality and polynomial computation time, these performance bounds are too broad, and far from the need of practical robotic deployment. Meantime, to obtain approximate results, some unrealistic assumptions or limitations have to be made. As a result, we conduct a comparative study on constructive heuristics, and find out the most promising heuristic that can produce good solutions in reasonable time.
In different heuristics, coupledness are dealt differently. In decoupled heuristics, two types of robots are separately scheduled. In coupled append or insertion, tasks are scheduled iteratively one by one. As of matching, multiple tasks are scheduled iteratively. Except the coupled insertion heuristic, other heuristics do not change the

Decoupled, greedy heuristic
For the multi-robot system involving two types of robots, decoupled or hierarchical planning are commonly used in most previous robotics literature [13,33,34]. Decoupled methods mean that firstly making decisions for one robot team, and then for the other robot team, with already decided variables acting as constraints. Using this idea, decoupled, greedy heuristic is designed as follows. 1) For each transport robot, find its subtask execution order separately by solving a traveling salesman problem, with the minimization of its routing component. This could be done by using existing traveling salesman problem solvers, either heuristics or metaheuristics. 2) Given all transport robots' subtask orders, pick robots' task scheduling problem can be regarded as a precedence-constrained task scheduling problem [35] belonging to the cross-schedule category in iTax. Pick robots' subtask orders are constructed iteratively by three rules, without violation of previously established transport robots' subtask orders.
(a) Free subtask first rule. Only subtasks that have no precedent subtasks in transport robot orders (i.e., "free" subtasks in [35]) are released. (b) Busy rule. Pick robots are not allowed to be idle if they have unfinished free subtasks. (c) Shortest rendezvous time first rule. When more than one free subtasks are in the work zone, the pick robot always chooses the subtask that requires minimal rendezvous time for the transport-pick robot pair. Free subtasks that have been appended to the pick robot subtask orders are removed form unscheduled subtask set.
The free-select-remove process continues until all subtasks are scheduled. Formal Analysis: We use an ordinary matrix in Eq.
(2) to analyze the time complexity of all algorithms. The dimension of an ordinary matrix transformed from a block rank matrix of Problem 1 is m × nK max , in which K max = max 1≤j≤n max 1≤i≤m K ij is the maximum block width. Then in the worst case each transport robot has at most nK max subtasks, and each pick robot has at most mK max subtasks.

Proposition 1 The decoupled greedy heuristic terminates in polynomial time.
Proof In step (a), suppose the traveling salesman problem is solved by the famous Christofides heuristic, which is an O (nK max ) 3 and 1.5-optimal algorithm. Totally the time complexity in step 1 is O mn 3 K 3 max . To generate a complete schedule, the free-select-remove process is run for |T | ≤ m × nK max iterations. Note that here |T | denotes the number of subtasks, rather than the number of tasks. In each iteration, at most m subtasks are released for n pick robots. The total time complexity for step (b)(c) is simply at most O(mnK max · mn). Finally, the worstcase time complexity of the decoupled greedy heuristic is O mn 3 K 3 max + m 2 n 2 K max . Then this proposition is proved.

Coupled, tripartite matching heuristic
The matching heuristic could generate rank-minimal task schedules. These rank-minimal schedules are of great interest because they could generate good task schedule in which the number of robots in the longest robotsubtask chain is minimized. The coupled tripartite matching heuristic works as follows. Let G * = (R T ∪R P ∪T , E * ) be the complete tripartite graph with the set of transport robots, pick robots and subtasks as three partitioned vertices. Each edge e ∈ E * is assigned with a weight that equals to the operate time of the subtask, plus the travel time between the robot's current vertex and the subtask vertex. In the tripartite matching, we consecutively determine a maximal number of subtasks (i.e., transportpick robot pair and mutual subtask) having the same rank (increasingly from 1 to r max ). In each iteration, at most max{m, n} subtasks are determined. This is done by solving the weighted tripartite maximum cardinality matching problem. Using this rule, the rank-minimal schedule is generated, i.e., the greatest value of the rank of the sequence is minimal.
In each iteration, based on how to determine the tripartite matching, different objectives may be used and thus many variants exist. We use the dynamic partial makespan as objective function, as Eq. (3) shows. The main idea of MINDYMAX is to choose robots and subtasks with minimal partial makespan in each iteration with all scheduled and chosen subtasks considered.
MINDYMAX : min max i,j,k∈M ha ijk + τ ijk , i = 1, 2, · · · , m, j = 1, 2, · · · n, k = 1, 2, · · · K max (3) Formal Analysis: It is well-known that solving a tripartite matching problem is generally NP-hard, so this heuristic may run infinitely based on the problem structure. However, as we do not need to solve the tripartite matching on complete tripartite graphs, the matching can be solved optimally in reasonable time. Meanwhile, as more and more subtasks are scheduled, the size of the tripartite matching problem decreases. As is shown in the latter empirical study, the coupled tripartite matching could generate task schedules in half minute.

Decoupled, bipartite matching heuristic
The decoupled, bipartite matching heuristic combines the main idea of two previous heuristics: making decisions for two types of robots separately; and obtaining a rank-minimal task schedule by matching. The first step is the same as step (a) in the decoupled, greedy heuristic. Because transport robots' subtask execution orders are pre-determined and do not change during scheduling, in each iteration we solve a bipartite maximum cardinality matching problem and at most max{m, n} subtasks are determined.
Formal Analysis: The algorithmic properties of the decoupled bipartite matching heuristic is summarized as follows.

Proposition 2 The decoupled, bipartite matching heuristic terminates in polynomial time.
Proof In the first step, the time complexity is O((nK max ) 3 ) as is shown in the proof of Proposition 1. The decoupled, bipartite matching heuristic terminates in at most 2d − 1 iterations [36], in which d = max{mK max , n}. The matching process is executed at most 2d − 1 times. The bipartite matching can be determined by the famous Hungarian algorithm, and its time complexity is O((max{m, n} · K max ) 3 ). Then the overall worst-case time complexity is O(n 3 K 3 max + 2 max{mK max , n} · (max{m, n} · K max ) 3 ). The proof is completed.

Coupled, append heuristic
In the coupled append heuristic, new subtasks are selected, scheduled and appended to the partial task schedule iteratively. Starting from an empty schedule, subtasks with the minimal start time of all unscheduled subtasks are appended. When multiple subtasks are available, the following three priority dispatching rules can be used as tie-breakers: STD (shortest time difference between transport-pick robot pair), SPT (shortest operate/processing time), and LPT (longest operate/processing time). Among these rules, SPT and LPT are commonly seen in shop scheduling, and STD is designed specifically for the TPS problem. As the append heuristic is easy to realize, all these rules can be used simultaneously and the best schedule is chosen. For the partial schedule, the makespan objective is computed exactly by the maximum of completion times of all scheduled subtasks.
Formal Analysis: The following proposition summarizes the time complexity of the coupled, append heuristic.

Proposition 3 The coupled append heuristic terminates in polynomial time.
Proof To generate a complete schedule, the append operation is run for |T | ≤ m × nK max iterations. In the p-th iteration, one subtask should be chosen from |T − p| by a particular sorting algorithm. Generally time complexity of the sorting algorithm for |T − p| elements can be estimated as O((|T | − p)log(|T | − p)). To sum up, the total worst-case time complexity for the coupled, append heuristic is O(m 2 n 2 K 2 max log(mnK max )). Then this proposition is proved.

Coupled, insertion heuristic
In the coupled insertion heuristic, for the p-th iteration, subtasks are successively inserted into a partial sequence S p according to a certain order, with all previous precedence relations preserved. To determine which subtask to be inserted, an insertion order can be determined by priority dispatching rules, or simply using previous heuristics. The makespan objective of a partial schedule is estimated by the maximum of the completion times of scheduled subtasks. For subtask t ijk , let v i and u j be the numbers of scheduled subtasks of the transport robot R T i and the pick robot R P j respectively. To insert the subtask t ijk into S p , (v i + 1)(u j + 1) possibilities exist. However, not every insertion leads to a feasible block sequence graph, because the sequence property has to be satisfied. The following three cases with inserting operations can lead to feasible solutions: C1 Let s ijk = 1, i.e., choose t ijk as the first subtask for both transport robot R T i and pick robot R P j . C2 Let s ijk = b + 1, where entry b occurs in the i -th row block or in the j -th column block of block rank matrix S p , i.e., subtask t ijk is a direct successor of one of scheduled subtasks of transport robot R T i or pick robot R P j . C3 If there exist a pair of subtasks t ipq and t rjs such that there does not exist any path between them in the current partial block sequence graph: -If r ipq ≤ r rjs , set s ijk = r rjs + 1 and r ipq = r rjs + 2, i.e., subtask t ijk is chosen as the direct successor of subtask t rjs and as the direct predecessor of t ipq . -If r ipq ≥ r rjs , set s ijk = r ipq + 1 and r rjs = r ipq + 2, i.e., subtask t ijk is chosen as the direct successor of subtask t ipq and as the direct predecessor of t rjs . -If r ipq = r rjs , both two previous rules should be used.
In all three cases, the ranks of all successors of the inserted subtask should be updated. Formal Analysis: Cases C1-C3 are the only cases that can lead to feasible sequence graph when inserting operations. The correctness and completeness of the coupled insertion heuristic are stated as the following two propositions.

Proposition 4 Cases C1-C3 generate all feasible sequences with the new subtask added to the current partial sequence.
Proof The proof is similar to the proof of insertion operation for open shop problems in [22,37]. Let SR(i) and SC(j) be the set of subtasks in row block i and column block j respectively, i.e., SR(i) = {s ijk |s ijk ∈ S}, SC(j) = {s ijk |s ijk ∈ S}. If we have a sequence on partial subtask set S, then we can surly construct a sequence on set S ∪ {s ijk }. Altogether we have (1 + |SR(i)|)(1 + |SC(j)|) possibilities to include s ijk in the new sequence. Now we have to eliminate all cyclic cases. The number of these cases is |SR(i)||SC(j)|−|A 1 ∪A 2 |, which means the number of new feasible sequences is: For all three cases, the number of new constructed sequences is: With Eqs. (4) and (5) this proposition is proved.

Proposition 5 The coupled, insertion heuristic enumerates all feasible sequences completely.
Proof This proposition is obvious. From an empty task sequence, in each iteration, all feasible partial task sequences are preserved. When all subtasks are inserted, all feasible sequences are recorded.
The following proposition shows the time complexity of the coupled, insertion heuristic.

Proposition 6 The coupled, insertion heuristic terminates in polynomial time.
Proof The insertion operation is executed for |T | ≤ mnK max times. Now we analyze the time complexity when inserting subtask t ijk into a partial block rank matrix S p in the p-th iteration. To update all successors of a changed subtask's rank in S p+1 with maximal rank r max p+1 , breadthfirst search should be used to analyze all its successors, and its time complexity is O(V p + E p ), in which V p and E p are the vertex set and the edge set of the partial block sequence graph S p . To obtain feasible block rank matrices with subtask t ijk included in the partial block rank matrix, C1 and C2 generate (1 + v i + u j ) new block rank matrices, and the time complexity is O((1 + v i + u j )(V p + E p )). In C3, to obtain a list of disconnected subtasks, for each of (v i + u j ) subtasks, breadth-first search should be used to analyze both its successors and predecessors, and the total time complexity is O(2(v i + u j )(V p + E p )). At most (v i + u j ) 2 disconnected subtask pair can be obtained. for each of them, to update all successors of changed subtask's rank in S p+1 with maximal rank r max p+1 , the time complexity is O((v i + u j ) 3 (V p + E p )). Totally for C1 to C3, the worst-case time complexity is approximated as O((v i + u j ) 3 (V p + E p )). For at most (v i + 1)(u j + 1) new feasible sequences, to evaluate their partial makespan, the time complexity is O((v i +1)(u j +1)mnK max ). Finally in the p-th iteration, the worst-case time complexity is approxi- As the insertion operation proceeds, v i → mK max , u j → nK max , V p → mnK max , E p → m 2 n 2 K 2 max . Then the previous time complexity can be rewritten as O((m + n) 3 m 2 n 2 K 5 max ) with lower order polynomials neglected. Totally, in the worst case, the whole coupled insertion heuristic will cost O (m + n) 3 m 3 n 3 K 5 max time. The proof is completed.
The following conclusions can be drawn from previous theoretical analysis. For the coupled insertion heuristic, its computation time is highly dependent on the number of scheduled subtasks in the partial task schedule. In order to guarantee non-redundancy and feasibility of new partial sequences, complex pre-processing subroutines are required. Generation of new feasible sequences are the most time-consuming subroutine. As more and more subtasks are scheduled, the computation time for inserting new subtasks grows notably.

Performance evaluation
To evaluate the performance of these heuristics, we present computational results obtained from our MATLAB-based simulation platform. Algorithms are run on a 2.9 GHz AMD Ryzen 7-4800H PC with 16 GB RAM. The attached also video shows the simulation platform and how the heterogeneous robotic order fulfillment system works.

Dataset design
For TPS problems, a comprehensive, discriminative and balanced evaluation is extremely challenging: a complete enumeration of all problem constituents, including environment maps, robot kinematics, numbers of two robot types, start locations, subtask distributions, number of subtasks in one task, number of subtasks in each zone, operate times, routing-scheduling ratio, and zoning strategy is unrealistic. Pragmatically, we concentrate on the most concerned constituents, i.e., robot heterogeneity, multi-station order fulfillment tasks, and simplify other less concerned constituents. Table 2 shows instance types in our dataset. Generally, an instance can be described as m-n-p-α, i.e., transport robot number, pick robot number, subtask number, and temporal-spatial ratio. For each type, we randomly generate 50 instances, and 1800 instances in total. Generation of instances are explained as follows.

Environment
All instances are generated on the same warehouse environment as Fig. 3 shows. In this warehouse (dimension 52 m × 39 m), 36 racks (dimension 10 m × 1 m) are arranged into a 9 × 4 array, thus generating 8 × 4 = 32 array of aisles for pick robots. These aisles can be divided equally into n work zones and allocated to pick robots. The aisle width is set to 3 m, a minimal number to allow robots go through when two robots occupying opposite stations, yet not causing waste of space. The width of cross-aisles and hall area all are 2 m, a minimal number to allow robots to move in two directions. Pickup locations are adjacent to racks. Delivery locations are placed in the periphery. This map is adequate for us to generate miscellaneous instance types with different zoning strategies.

Robots
We assume a holonomic disc robot model with unit speed 1 m/s. As is shown in Fig. 4, to avoid collisions when one robot goes to a vertex where the other robot departs from, the diameter of the robot should be less than ≤ √ 2 2 m. All robot share the same kinematics, and they are discriminated only by their task roles. In each discrete time step, a robot performs an action in the action set {↑, ↓, ←, →, }, i.e., move up, down, left, right, and stay still. Transport robots randomly start from delivery stations, and all pick robots start from fixed initial stations. Transport robot numbers are among 5,10,15,20,25,30. Pick robot numbers are among 8,16,32, with zoning strategies of 2×2 aisles,1×2 aisles,1×1 aisle respectively.

Tasks
In one task, we fix the number of subtasks as 10, a moderate task size. Comparatively in open shop scheduling problems, number of operations in each job is the machine number. However, we cannot use this rule because realistic order fulfillment tasks are surely not related to pick robots. We also notice that random numbers (e.g., uniform distribution U (5,15)) are better, however, this would make it difficult to control instances, because we have allowed subtasks to be located randomly in any work zones. In other words, we stress randomness more on subtask spatial distributions rather than number of subtasks in one task. We have to make these compromises although the introduced heuristics can handle any circumstances. In one task, pickup subtasks are randomly located in pickup stations. Delivery subtasks are randomly located in delivery stations. For small (α = S) and large (α = L) temporal-spatial ratios, the operate time is generated subjecting to uniform distributions U (10,20) and U(10, 50) respectively. Each task contains a pre-optimized subtask execution order with minimized routing components obtained by an optimal traveling salesman problem solver. Note that for instances with different pick robot numbers, we generate the same tasks, so that we can explore relations between the overall system performance and different transport-pick robot number combinations. In other words, totally 600 tasks are generated.

Methods
The proposed constructive heuristics are abbreviated as follows: • DG : Decoupled heuristic.
• CI : Coupled, insertion heuristic, with insertion order generated by append heuristic.
In practice, actual travel time between two different subtask stations can be estimated by naive multiplicative robustification, i.e., multiply the ideal travel time by a factor larger than 1. Here we simply use the ideal travel time, obtained by the ideal geodesic distance between them with obstacles considered, divided by unit speed. This is feasible, because the ideal task schedule can always be scaled by the same factor and transformed to the task schedule when multiplicative robustifucation is used. For each heuristic on each instance, we record its normalized performance ratio and computation time. The normalized performance ratio obtained by a certain heuristic * is computed by C max ( * )/LB.

Results and analysis
First of all, for append heuristic and different priority dispatching rules, we count how many times a rule generate best value, as Table 3 shows. For instances with n = 8, STD works best; in contrast for n = 16, 32, LPT works best. This is because for smaller pick robot numbers, both pick robots' and transport robots' routing components are dominant, and STD rule could reduce the time difference that is required for the transport-pick robot pair rendezvous. For larger pick robot numbers, pick robots' work zones are getting smaller, and their routing component are less influential. In latter result analysis, we only choose append heuristic whose priority dispatching rule generate the best value as a representative. For n = 8, we choose STD, and for n = 16, 32, we choose LPT. Figure 5 shows the normalized performance ratio on different instance types. First, DG and DM generate similar results with largest performance ratio, largely because of their restrictions on problem search space. CM performs better than DG and DM, except for instances with n = 32, p = 50, 100, 150 in m-32-p-S instance types and n = 32, p = 50, 100 in m-32-p-L instance types. This is because with few subtasks, few transport robots and more pick robots, transport robots' routing components dominate the overall optimization process. In DG and DM these components are firstly minimized. For all other cases, routing and scheduling components of both robot types contribute equally to the overall optimization. The normalized performance ratio for CA and CI closely intertwine and outperform all the other heuristics regardless of problem types.

Optimality
Specifically, the performance gap between coupled CM, CA, CI and decoupled DG, DM are more obvious on instances with relatively few pick robots (n = 8, 16), as pick robots' work zones are larger, and its routing subproblems are more influential. As more pick robots are used, work zones become smaller, and their routing components have a less influence on the whole problem. Extremely speaking, there is a tendency for pick robots to be stationary rather than mobile.

Computation time
We only report computation times for instances with α = S (results for instances with α = L are similar), as Table 4 shows. Computation times for DG, DM, and CA increase polynomially with lower-degree. For CI, computation times increase tremendously. It is observed that in the insertion process, as more subtasks are scheduled, more computation time is also required. For DG, DM, CM, and CA, with large pick robot numbers, more computation time is required. Oppositely for CI, with small pick robot numbers, more computation time is required. This is because in these instances, pick robots' work zones are larger and they have to fulfill more subtasks. The number of scheduled subtasks of pick robots grows more rapidly than the other instances. Although both append and insertion heuristics generate similar results, append is more desirable because it is simple to use and fast to compute.

Overall system performance and robot numbers
As we generate the same tasks for different pick robot numbers, we can analysis the relationship between the overall system performance and robot numbers. Take optimality results obtained by append heuristic for example. Adding more pick robots (from n = 8 to n = 16 and from n = 16 to n = 32) will definitely improve the system performance, however, the averaged improvement per pick robot will become less obvious as pick robot number increases. This is subject to the law of diminishing marginal utility [38] in economics. Practitioners can use this phenomenon to balance balance total cost (i.e., numbers of two types of robots) and throughput requirement.
Remark We notice that the computational results and performance analysis in this section are surely not comprehensive enough, because many other problem constituents are neglected. This is a compromise that we have to make, because it is impossible to cover all combinations. Meanwhile, currently there is also a lack of benchmarks. Nevertheless, the coupled append heuristic can  always be recommended for all circumstances, because of its superior optimality, and lower time complexity order. It is simple to realize, and previous partial task schedules are unchanged.

Conclusions and future work
In this work, we conduct a comparative study on constructive heuristics using different principles for the TPS problem. Theoretically we provide formal analysis of these heuristics. Empirically, we draw the following conclusions from computational results. 1) Append heuristic is highly recommended for the TPS problem, with better performance in most cases, and within reasonable computation time. 2) Specifically, coupled heuristics work better on instances with relatively few pick robots and large work zones. 3) On the overall system performance and robot numbers, the law of diminishing marginal utility is observed.
There are many possibilities for future work. On task scheduling, we merely present several fundamental heuristics, in fact, more combined heuristics and metaheuristics methods could also be used. Our research methodology is essentially empirical, and there is urgent requirements on developing polynomial approximate methods with bounded suboptimality. More realistic constraints could also be considered, such as task release times and due times, robots' energy constraints, robot failure. Meanwhile, online scheduling with unknown operate times and uncertain travel times is also paramount. On task execution, to enable effective and robust execution of planned task schedules, developing simultaneous task and path planning algorithms is necessary. On heterogeneous robotic order fulfillment systems, future work may include decision-making in strategic level, tactic level and operational level, including general principles on system and layout design, numbers of two types of robots, long-term system performance with time-changing workload.