Scheduling a pick and place packaging line

In this paper, we introduce the Pick and Place Packaging Problem (P4), for optimally scheduling a packaging system with one input conveyor (pick) and one output conveyor (place). We give a formal definition of the underlying dynamic optimization problem. We describe two properties that hold for its solutions and present an efficient row generation approach exploiting these properties. Starting from a basic version of the problem, we introduce two variants where we account for the possibility of holding the grip of items and variable conveyor speed. We extend the proposed exact solution method to these two cases. Then, we present an efficient Iterated Local Search heuristic for the problem and its variants. Numerical results show the effectiveness of our approach.


Introduction and state of the art
We consider the problem of scheduling the operations of a Pick and Place system in a Packaging line. Pick and Place systems consist of a set of robotic arms and one or more conveyors. One of the conveyors serves as the input of the system; thus, we refer to it as the "input conveyor." The robots are tasked to remove the items flowing on the input conveyor and place them in a designated destination. All items that reach the end of the conveyor without being picked up are lost. A dynamic example of this type of systems can be watched on the site: https://www.youtube.com/watch?v=PrZw-xp_VzE In this paper, we consider the case of a Pick and Place line with two conveyors, one input conveyor and one output conveyor. The robots have to place the input items in some B containers in the output conveyor. The items are identical, and the containers have a limited capacity. The operation of picking up one item from the input conveyor and placing it to the final destination is called mission. The problem that we study here involves making decisions about the scheduling of the missions and, in the extension of the problem, also about the conveyor's speed. The objective function consists of the minimization of the number of items that are not picked up. In such problems, efficiency is of paramount importance. Utilizing the robots to their fullest potential is essential due to their high cost.
We assume that the position of the incoming items is known in advance. Generally, in these systems, highdefinition cameras monitor the inbound conveyors shortly before the items enter within the reach of the robots and provide information of the location of all incoming items. Typically, these cameras cover a few meters before the reach of the first robot. The position of the input items is usually known a very short time in advance (i.e., 5-10 seconds). Thus, the scheduling of the robots has to be decided in a timely manner. For this reason, the majority of methods proposed for this kind of system implement simple dispatching strategies, and the most commonly used approaches in the industry use queuing mechanisms. These techniques perform rather well when the system reaches the steady-state conditions. However, in systems with frequent transients and irregularities in their operating conditions, this kind of strategy often falls short, and the performance gets signifi-cantly deteriorated. The alternative to these strategies is to repeatedly re-optimize the scheduling of the line as more information becomes available. This is generally referred to as a rolling horizon approach, where the scheduling periodically is recomputed for the next planning horizon. Then, when the computed scheduling is applied, the scheduling of the next period is solved.
Most papers on the optimization of Pick and Place systems concern themselves with the optimization of single conveyor systems and propose techniques that do not generalize well in the two-conveyor case. In these cases, the optimization effort considers maximizing the system throughput and minimizing the wear. This comes in two main forms: -The load balancing of the various machines, so that no machine is ever overworked; -The minimization of the robot's travel time, by optimizing of the mission starting time.
In their papers, Mattone et al. (1998) and Mattone et al. (2000) presented augmented versions of the commonly used FIFO and LIFO dispatching rules for a single conveyor Pick and Place system. A more advanced approach is proposed by Huang et al. (2012), it combines heuristic dispatching rules, to obtain a performant and fast approach. Humbert et al. (2015) study the performance of several approaches mostly focusing on the various developed dispatching rules in an experimental setting. Huang et al. (2015) present part dispatching rules in a task-based environment. Bozma and Kalalıoglu (2012) propose a coordination approach based on cost-driven policies. These advanced methods allowed to directly pursue more elaborate objectives, and the heuristic and metaheuristic approaches ensure fast computations. For the minimization of robot travel times on a simple conveyor system (Daoud et al. 2010) concentrated on the computation speed using a metaheuristic approach. Ahmadi and Kouvelis (1994) put their attention on a different aspect of the problem, trying to balance the usage of the different robots on the line as to minimize their wear. The energy optimization and the minimization of the damage of the components are studied by Pellicciari et al. (2013). Ho and Ji (2009) focus on the stationary problem that is solved exactly to minimize the traveled distance of the robots, as in a sequence-dependent scheduling problem.
The problem of scheduling two conveyors however, to the best of our knowledge, has received little interest over the years. Despite this, the problem poses interesting challenges, both on the optimization of the missions of the robots themselves, and on the synchronization of the various robots acting on the system. For the problem of controlling a packaging line for perishable Ferrari et al. (2015) and Pizzi et al. (2016) proposed the idea of using a rolling horizon approach paired with an exact model to dynamically compute the scheduling of the robots and the conveyor's speed. They were able to show significant benefits for the control of the line over the commonly used dispatching strategies. This approach is very promising, though it can be applied only to very small instances, due to the cumbersome scheduling computation.
In this paper, we give a formal definition of the scheduling problem proposed in Ferrari et al. (2015). We will refer to this problem as the Pick and Place Packaging Problem (P4). We will present an exact solution method for a basic version of the problem. This approach will be extended to consider the case when the robot can hold the grip of items, and the case when the conveyor speed can be varied. In order to cope with real instances in a manageable time that is compatible with a rolling horizon approach, we will propose an efficient heuristic based on an Iterated Local Search framework. We will test the developed algorithms on a set of instances based on three realistic lines to validate the developed solution method and Iterated Local Search heuristic.
The rest of the paper is organized as follows: in Sect. 2, we give the formal description of P4 and define two theoretical properties of the problem that can be exploited in its solution. We propose an efficient solution method based on row generation in Sect. 3. Then, we deal with the extensions of the problem to the case where robots can hold an item before moving it to the containers (Sect. 4.1) and to the case where the speed of one of the conveyors is controllable (Sect. 4.2). Afterward, we develop a heuristic method based on Iterated Local Search in Sect. 5. In Sect. 6, we present the conducted computational experiments on the problem. Lastly, in Sect. 7 we draw our conclusions and discuss future research.

Problem details
The Pick and Place Packaging Problem (P4) considers a system with several robots and two conveyors, denoted as the " Input " and the " Output " conveyors. In Fig. 1, we show a schematic view of one possible line. The two conveyors are parallel and move in the same direction. However, the methods presented here can be extended to systems with any configuration since they not require specifically this assumption. The velocity of the conveyors is fixed. In Sect. 4.2, we relax this assumption and allow our method to control the velocity profile of one of the conveyors.
The items to be handled by the system enter the line through the Input conveyor. The task of the system is to pick up those items and move them in one of the containers on the output conveyor. All items are identical, aside from their position on the input conveyor. We consider the items partitioned into columns, i.e., groups of items occupying the same longitudinal position on the conveyor. Let I be the set of the . . , n}. We denote by n i the number of items in column i at the beginning of the planning horizon before any item is picked by the robots.
Let C be the set of the containers on the output conveyor C = {1, . . . , m}. We denote by m j the maximum capacity of container j at the start of the planning horizon. For the sake of brevity, we will use " item i " as a short-hand to refer to an item of column i, with i ∈ I , and " slot j " as a short-hand to refer to one of the empty spots in container j ∈ C.
The operation of moving the items to their destination is handled by the robots that are present along the line. We denote as R the set of robots. The robots are identical, singlegrip robots, with non-overlapping working areas. Thus, the robots can move one item at a time, and the operations of each robot do not interfere with the other robots on a mechanical level.
A single pick-and-place operation performed by the robots is called " mission ". We use the tuple (i, j, r ) to identify the mission where robot r ∈ R picks up an item from column i ∈ I and puts it into container j ∈ C.
A solution of the P4 involves three levels of decisions: (1) determining the set of missions performed by the robots, (2) determining the order in which the missions are performed, and (3) assigning a starting time to each of those missions.
Let us define an " assignment " A as the set of missions of the solution, and A r ⊆ A as the assignment of robot r , i.e., the set of missions assigned to r in A. We denote as A r i, j the number of missions (i, j, r ) in A, ∀i ∈ I , ∀ j ∈ C, ∀r ∈ R. An assignment is valid, if: Given A r , we define a " sequencing " S r A as the ordered sets of missions performed by each robot r . We denote by S r A [ p] the p-th mission of S r A . Additionally, given A, we define a " scheduling " T A as a set of starting times t (i, j,r ) of each mission (i, j, r ) ∈ A, and let T r A denote the scheduling of A r . We note that a scheduling of a given assignment implicitly entails a sequencing for the same assignment. When a sequencing S r A corresponds to the sequencing entailed in a scheduling T r A , we say that that T r A "matches" S r A .
When A contains multiple copies of the same mission (i, j, r ) (i.e., A r i, j ≥ 2), (i, j, r ) p indicates the p-th mission (i, j, r ) in A. We denote by t p (i, j,r ) the starting time of mission (i, j, r ) p . The problem input is given by the following data, for each robot r , column i and container j: a r i and b r i are the earliest and latest time at which item i can be picked up by robot r , -α r j and β r j are the earliest and latest time at which one item can be placed in container j by robot r .
We assume that if a r i < a r i then b r i < b r i , and if α r j < α r j then β r j < β r j . Without loss of generality, we assume that the indices of I and C are in order of earliest starting time. Thus, ∀i, i ∈ I : i > i ⇒ b r i > b r i and a r i > a r i ∀r ∈ R, and ∀ j, j ∈ C : j > j ⇒ β r j > β r j and α r j > α r j ∀r ∈ R. A mission (i, j, r ) is feasible if its starting time falls within its availability window, i.e., (1) We initially assume that the time for performing a mission is constant and requires δ units of time. This implies that considering mission (i, j, r ), if there is no intersection between the availability windows of item i and slot j for robot r , the mission is infeasible. However, in some systems, the robot may pick up an item and hold it until the destination slot becomes reachable. This would extend the set of possible feasible solutions but will introduce a variable mission time δ (i, j,r ) . For the sake of simplicity, we discuss our exact solution method for the simple case with constant mission time to extend it afterward to the most general case.
Let E r denote the set of the possible missions of robot r ∈ R, i.e., the set of missions (i, j, r ) A is feasible if all the missions in A r are feasible, and the starting times any two distinct missions (i, j, r ) p , (i , j , r ) p ∈ A r are separated by at least the time needed to perform the first mission, i.e., If there exists a feasible scheduling T r A for assignment A r , we say that the assignment A r is feasible. Similarly, if there exists a feasible scheduling T r A that matches a sequencing S r A , we say that S r A is feasible. Lastly, if an assignment A is feasible for all robots, we say that it is feasible.
If an item is not picked up by any mission in the assignment, the item is lost. The number of lost items is given by: The P4 problem consists of finding an assignment A and a feasible scheduling T A , such that the total number of lost items is minimized. An alternative objective function is minimizing the number of lost slots, i.e., slots that are not filled by any mission in the assignment.

Properties of the solutions
The solution methods that we propose are based on two properties.
Property 1 Let A r be an assignment for robot r , such that Proof Assume without loss of generality that robot r performs the mission (i, j , r ) before (i , j, r ) (see Fig. 2). Additionally, for simplicity, assume that missions (i, j , r ) and (i , j, r ) are performed immediately one after the other, with no idle time.
Given that A r is feasible, it holds: On the other hand,Â r is feasible if: We can then rewrite the conditions that must hold for the feasibility of A r as: Similarly,Â r is feasible if: Recalling that b r i < b r i , a r i < a r i , β r j < β r j , and α r j < α r j ; conditions (3)-(4) imply conditions (5)-(6).
Definition 1 An assignment A r having no two missions (i, j , r ), (i , j, r ) with i < i and j < j , is" 1-irreducible ".

Corollary 1 Let A r be a 1-irreducible assignment for robot r, if A r is not feasible then the assignment:
Property 2 Let A r be a 1-irreducible assignment for robot r , such that (i, j, r ), (i , j , r ) ∈ A r , with i ≤ i and j ≤ j , and let there be a feasible sequencing S r A such that (i , j , r ) precedes (i, j, r ). Then, there exists another feasible sequencinĝ S r A of A r such that (i, j, r ) precedes (i , j , r ). Proof Let us denote by T r A a feasible scheduling that matches S r A . Additionally, suppose thatT r A is a scheduling that matchesŜ r A where missions are performed as early as possible. We now show that if T r A is feasible,T r A is also feasible. Without loss of generality, let us assume that in S r A andŜ r A all missions are performed as early as possible. Lastly, for simplicity, let us assume that robot r is idle before performing missions (i, j, r ) and (i , j , r ) (see Fig. 2). Under these assumptions, we can write the starting time of the two missions (i , j , r ) and (i, j, r ) in the two schedulings. In case T r A , it holds t (i , j ,r ) = max(a r i , α r j ), and t (i, j,r ) = max(t (i , j ,r ) + δ, a r i , α r j ). We explicitly write the feasibility conditions of T r A as such: For schedulingT r A , it holds t (i, j,r ) = max(a r i , α r j ) and t (i , j ,r ) Conditions (7)-(8) and (10)-(11) coincide. Conversely, condition (9) implies (12). Thus, if S r A is feasible,Ŝ r A is also feasible.

Definition 2 A sequencing S r
A is "2-irreducible" if there are no two missions (i, j, r ), (i , j , r ) ∈ A r , with i ≤ i and j ≤ j , at least one of which holds strictly, such that (i, j, r ) follows (i , j , r ).

An exact method based on a row generation framework
We now present a MILP formulation for the P4. We express using variable x r i, j ∈ Z + the number of times mission (i, j, r ) is performed in the scheduling, i.e., x r i, j = A r i, j . Binary variable y r i, j is one if x r i, j > 0; and to zero otherwise. Variable r i, j is the starting time of mission the first mission (i, j, r ) performed in the schedule, i.e., = r i, j = t 1 (i, j,r ) . If no mission (i, j, r ) is performed (i.e., y r i, j = 0), variable r i, j is free to take any value. As for Corollary 2, without loss of generality, we assume that when multiple missions (i, j, r ) are scheduled, those missions are performed consecutively starting at time r i, j . Thus, t p (i, j,r ) = r i, j + ( p − 1)δ. Lastly, variable λ i is the number of items in column i ∈ I that are not assigned to any mission (i.e., the number of items lost). The formulation is as follows.
The objective function (13) minimizes the number of items lost. Constraints (14)-(15) ensure the feasibility of all scheduled missions with respect to the availability of items and containers. Constraints (16) guarantee that the assignment is 1-irreducible. Constraints (17) prevent the robots from performing more than one mission at a time. In those constraints, M is an arbitrarily large constant. Constraints (18) and (19) enforce the availability of items and the capacity of the containers. Constraints (20) link the y r i, j and x r i, j variables. Lastly, constraints (21) and (22) are the domain constraints. Note that if no mission (i, j, r ) is performed (i.e., y r i, j = 0), the value of r i, j has no significance to the solution. In this case, constraints (14)-(15), and (17) can always be satisfied by setting r i, j = max(a r i , α r j ). We develop a dynamic constraint generation algorithm to solve the formulation. Specifically, our algorithm can be summarized in the following steps: 1. Generate and solve a Master Problem (MP) to obtain an initial assignment A. 2. Apply property 1 to the initial assignment to obtain an equivalent 1-irreducible assignmentÂ and compute a 2irreducible sequencing SÂ. 3. Generate TÂ and verify its feasibility. 4. If TÂ is feasible, returnÂ and TÂ. Otherwise, separate one or more violated cut for the MP and iterate the procedure.
All the steps will now be described in detail.

Master problem
The Master Problem is used to generate an initial assignment A, accounting for the availability of items and containers. The MP is formulated as follows.

Refinement step
In the refinement step, we apply property 1 to A, to obtain an equivalent 1-irreducible assignmentÂ and compute a 2irreducible sequencing SÂ. This procedure produces a set of missions that dominates any other one when attempting to generate a scheduling of the assignment as for Corollary 1 and Corollary 2. Moreover, a 1-irreducible assignment is needed to generate tight constraints in the next phase. The refinement step is carried out separately for each robot. For a given robot r ∈ R, we store in vector I r the number of items of each column i ∈ I assigned to r ∈ R in A. Similarly, we use vector C r to store the number of slots of each container j ∈ C assigned to r ∈ R in A. We con-structÂ r and S rÂ iteratively from an empty assignment and We append (i , j , r ) to S rÂ and decrease both I r [i ] and C r [ j ] by one. This process is iterated until all elements of I r and C r are zero. This procedure is described by Algorithm 1, whose time complexity is linear in the number of missions.
Algorithm 1: The procedure used for refining the assign-

Time feasibility check
After the refinement step of each robot r , the feasibility of S rÂ is checked in terms of time. The feasibility check is carried out separately for each robot, by considering the missions of the sequencing in the order, and iteratively assigning them a starting time. The process continues until either all missions have been considered, or an unfeasibility has been found. We assume without loss of generality that all missions are scheduled at their earliest possible time. Therefore (i, j, r ) will start as soon as the robot r is idle after the previous mission, and the column i and the container j are both available to r . Let us denote by ψ p (i, j,r ) the time instant when robot r becomes available for mission (i, j, r ) p . The ready time of where (i , j , r ) p is the mission immediately proceeding (i, j, r ) p in the scheduling of r . The starting time of (i, j, r ) p is obtained as follows.
If the choice of t p (i, j,r ) given by (27) is not compatible with (1), then there exists no feasible scheduling forÂ r as a consequence of corollary 2, nor for A r as for corollary 1. When this occurs, we terminate the time feasibility check for r .
The time feasibility check of a given robot is described by Algorithm 2, whose complexity is linear in the size of A r .

Algorithm 2: Time feasibility check for robot
If Algorithm 2 terminates having assigned a starting time to each mission,Â r is feasible. Conversely, ifÂ r is infeasible we generate a new constraint for robot r . Therefore, at each iteration, we can generate at most one constraint for each robot.

Dynamically generated constraints
When we detect an infeasible mission, we amend the MP generating an additional constraint. In the rest of this section, we first present a class of valid inequalities for the problem. Then, given an infeasible solution, we discuss how the parameters of the presented inequality are selected to generate a violated cut. Lastly, we present a lifting of the generated inequalities that produces stronger cuts for the problem.
For the sake of brevity, let us define a "section" as a subset of the possible missions of a robot contained between a starting and an ending missions ( f , F, r ) and (l, L, r ) that is: The latest feasible starting time of any mission in the section is t K = min(b r l , β r L ). Thus, the starting time of all the missions inĒ r has to fall within the time interval [t 1 , t K ). The maximum number of missions that can be contained inĒ r is given by: The following constraint constitute a valid inequality for assignment A: We now discuss the choice ofĒ r that yields a violated constraint upon detecting an unfeasibility inÂ r . Let mission (i , j , r ) be the infeasible mission detected inÂ r . Additionally, let (i , j , r ) be the last mission before (i , j , r ) inÂ r such that: According to (27), all missions in E r (i , j , i , j ) are performed with no idle time, and the starting time of (i , j , r ) can be computed as: Since t (i , j ,r ) violates (1), it results that: Consequently, settingĒ r = E r (i , j , i , j ) yields a violated constraint forÂ r . Writing (29) utilizing the MP variables leads to the following valid inequality: Despite being valid, constraint (30) is rendered ineffective by the refinement step that maps A r toÂ r . Indeed, it is possible to find an assignment A r that satisfies constraint (30), but that, refined toÂ r , violates it. To understand how this can happen, consider the following example.
We strengthened (30) by enlarging sectionĒ r toẼ r = E r ( f , l, F, L). Specifically, we resolve the min and max operators in (28) to determine the largest sectionẼ r ⊇Ē r , The resulting constraint added to the MP is the following: In our previous example, this would correspond to extending the section of the constraint from E r (1, 1, 2, 2) to either E r (1, 1, 3, 2) or E r (1, 1, 2, 3), depending the availability intervals of the items and boxes. Respectively giving rise to either constraint: or constraint: x r 11 + x r 12 + x r 21 + x r 22 + x r 13 + x r 23 ≤ 1.
Both of these constraints would prevent the cycling behavior presented. These " extended " parameters strengthen the constraint, reducing the number of times the MP needs to be solved and preventing the MP to bypass some of the added constraints by exploiting the refinement step.

Problem extensions
So far, we have assumed a simple problem setting. Namely, we assumed a constant mission time, deriving from the fact that the robots place the items immediately after they have picked them up; and a constant conveyor speed. This setting is enough representative when the input and output flows are reasonably steady. Indeed, if items and containers arrive on the conveyors with good regularity, placing items just after having picked them up is by far more efficient. However, when the items and the containers start to be highly scattered, the impossibility of holding the items by the robots to wait for the arrival of a container, or to modify the speed of one of the conveyors may be limiting. Thus, excluding some missions that may improve the quality of the scheduling. In this section, we propose two extensions. The first one allows the robot to hold the grip of the items. The second one allows us to modify the speed of one of the two conveyors.

Item holding
When a robot is allowed to hold an item, the assumption of having a constant mission time does no longer hold true. Indeed, a mission (i, j, r ) may start when the items i is available to r but the container j is not yet in the reach of the robot. Thus, formally, t (i, j,r ) is not required to be in the intersection of the two intervals [a r i , b r i ) and [α r j , β r j ). The feasibility condition (1) is extended as follows.
If i and j are such that a r i < α r j and t (i, j,r ) ∈ [a r i , α r j ) the mission is feasible; however, its duration will be larger than δ. The robot has to wait α r j − t (i, j,r ) before the j is reachable; thus, the total mission time is δ + α r j − t (i, j,r ) . We note that in this case, the mission finishes at time δ + α r j , regardless of the choice of t (i, j,r ) .
Without loss of generality, we assume that the starting time of the missions will be selected such that its duration is minimized. Thus, we distinguish two cases, either (1) have an empty intersection. In the first case, mission (i, j, r ) has a duration of δ when t (i, j,r ) satisfies (1). If the availability of i and j has an empty intersection, and β r j < a r i , mission (i, j, r ) cannot be performed. If b r i < α r j , the mission is feasible but will involve some waiting. The minimum duration of the mission is δ + α r j − b r i , and is achieved when t (i, j,r ) = b r i . Let us define δ (i, j,r ) as the minimum duration of mission (i, j, r ), i.e., Using δ (i, j,r ) , condition (2) is updated as follows.
It can be shown that properties 1 and 2 also hold in this extended case.

Extension of the formulation
In order to adapt the formulation to the case where the robots may hold an item during a mission, we explicitly distinguish the case where δ (i, j,r ) = δ (i.e., the holding case) from the case where δ (i, j,r ) > δ. We replace constraints (14)- (15) with the following ones: Constraints (38) enforce the feasibility of t 1 (i, j,r ) . Constraints (39)-(40) enforce the feasibility of all missions (i, j, r ) when δ (i, j,r ) = δ. We observe that δ (i, j,r ) > δ implies that in a feasible scheduling x r i, j ≤ 1. Constraints (41)-(42) enforce x r i, j ≤ 1 when δ (i, j,r ) > δ and fix the starting time of those missions to b r i . We update constraint (17) to consider the variable mission duration.

Extension of feasibility check
We adjust the time feasibility check to account for item holding by updating the computation of the starting time of the missions (27) as follows.
Additionally, the computation of the ready time of the robots is updated to use δ (i, j,r ) in place of δ.

Extension of the generated constraints
When the item holding is allowed, the earliest starting time of the first mission of a section is only bounded by the availability of the items, t 1 = a l . However, any mission that starts before the first container is available will have to wait. Thus, the earliest starting time of the second mission in the section is t 2 = max(a f , α F ) + δ, and the earliest starting time of the p-th mission is t p = max(a r l , α r F ) + ( p − 1)δ. In general, if (min(b r l , β r L ) − max(a r f , α r F )) > 0, the capacity of a section E r ( f , l, L, F) follows the same formula as the non-holding case (28). On the other hand, if (min(b r l , β r L ) − max(a r f , α r F )) ≤ 0, we need to distinguish two cases, β r L − a r f ≤ 0, and b r l − α r F <= 0. In the first case, the first item of the section becomes available after the last container is no longer within reach of the robot. Therefore, the capacity of the section is u r ( f , l, L, F) = 0. In the second case, the robot may pick up an item and wait for the first container of the section to become available; thus, the capacity of the section is u r ( f , l, L, F) = 1.
Lastly, we discuss the choice of the section on which to impose the generated constraints. We recall that when generating constraints, we consider the sectionĒ r = E r (i , j , i , j ), and later expand it using (31)-(34). Where mission (i , j , r ) is the infeasible mission detected in the current solution, and (i , j , r ) is the last mission before (i , j , r ) inÂ r such that: So that all missions in E r (i , j , i , j ) are performed with no idle time.
We update the choice of mission (i , j , r ) to account for the possibility of holding the items, as the last mission before (i , j , r ) ∈Â r such that:

Controllable output conveyor
We now consider the possibility of varying the output conveyor speed. We adapt the exact algorithm of Sect. 3 to this case.

Representation of the conveyor speeds
To represent the speed dimension, we apply a discretization of the conveyor length into fixed steps. Let {1, . . . , S} be the set of steps of length l dividing the conveyor space. We denote by τ s the time the conveyor takes to move along step s ∈ {1, . . . , S}, we refer to this quantity as " step time ". During that time, we assume that the conveyor maintains a constant speed given by v s = l τ s . See an example of discretization and speed profile in Fig. 5.
Using this representation, we can rewrite any distance d on the conveyor as d = γ l + σ l, with γ ∈ {0, . . . , S − 1} and σ ∈ [0, 1). Meaning that to cover distance d, the conveyor needs to make γ steps of length l and then cover a portion σ Fig. 5 An example of the resulting motion that can be obtained decoding of the speed profile. On the left the position of the conveyor over time, on the right the speed of the conveyor of the next step. Consequently, the time needed to cover d is: Using the position of the containers on the line, it is possible to rewrite the availability interval of a container as a function of the conveyor speed. Let the quantities γ r j , γ r j , γ r j , and γ r j describe the position of the containers on the conveyor. We express α r j and β r i as follows.

Extension of the method
To adapt the model and the solution method to the variable conveyor speed case, we redefine the parameters α and β depending on variables τ s using (45) and (46). For each generated constraint, we introduce an additional binary helper variable z, which equals one if the constraint is deactivated and zero otherwise. The generated constraints are modified as follows.
Where M is an arbitrarily large constant. To generate the modified constraints, we follow the same procedure described in the fixed speed case. Letβ r j andα r j be equal to the values that β r j and α r j assume under the optimal conveyor speed at the moment of the constraint generation. The parameters that define the constraint are set as follows: These constraints have been modified to handle three issues. Firstly, the section E r ( f , l, F, L) of constraints (35) is dependent on α and β. We use (50-53) to resolve the outcomes of (31-34) at the generation of the constraints.
Secondly, the expression of the capacity of a section (28) is nonlinear in the availability interval of the containers. We use (49) to resolve the outcome of the min and max operators in (28) beforehand.
Lastly, the capacity of a section can be a function of the availability of the containers. Therefore, the capacity of the constraints can become negative under some velocity profiles, which would make the problem infeasible. The constraint (47) can be deactivated to account for this occurrence and guarantee feasibility. In this case, the constraint (48) guarantees that no mission in the section is performed.
Finally, M is set to the maximum negative value that the capacity of the constraint can assume, that is:

An ILS based heuristic
The real-time nature of the problem and its complexity requires timely solution approaches that are not compatible with the exact ones. To cope with the requirements of a rolling horizon approach, we develop a metaheuristic algorithm for the problem based on an Iterated Local Search (ILS) approach (Lourenço et al. 2019). The advantage of using a metaheuristic is twofold. Firstly, metaheuristic algorithms are generally much faster than exact ones, making them generally more suited to real-time applications. Secondly, iterated heuristics can be terminated early and still provide a feasible solution, whereas this is not guaranteed to be the case for exact algorithms. The ILS alternates a local improvement phase and a perturbation phase to escape local optima. The algorithm terminates after N I L S iterations of perturbation and local improvement. The basic framework of Iterated Local Search is detailed in Algorithm 3.
return s best ; In the following sections, we will detail the components of the implemented ILS.

Solution representation
In the ILS, the solutions are represented as a separate allocation of items and slots to the robots. For each robot r ∈ R, vector I r stores the number of items allocated to robot r for each column i ∈ I . Vector C r stores the number of slots allocated to robot r for each container j ∈ C. To evaluate the solutions, we need to translate vectors I r and C r into a feasible assignment and scheduling. This decoding can be computed as the solution of the following optimization problem. Given the vector sets: I r and C r ∀r ∈ R, find the assignment A and scheduling T A that minimizes the number of lost items and that satisfies the following condition.
The decoding is carried out separately for each robot r constructing A r and T r A iteratively. At each iteration, we select mission (i , j , r ) as a candidate to be added to the assignment, where i = min{i : I r i > 0} and j = min{ j : C r j > 0}. We attempt to generate a starting time for mission t (i , j ,r ) following (27). If the generated starting time is compatible with (1), we add mission (i , j , r ) to A and t (i , j ,t) to T r A and decrease both I r [i ] and C r [ j ] by one. Otherwise, the candidate mission is infeasible. When this is the case, either the current column or the current container has to be discarded and considered lost. If b r i > β r j , we discard the current container j and we set C r [ j ] ← 0. Otherwise we discard the current column i and set I r [i ] ← 0. The process is iterated as long as both I r and C r have at least one nonzero element.
The pseudocode of the decoding procedure is presented in Algorithm 4.

Algorithm 4:
The procedure to decode the solution of

Local improvement phase
The local improvement phase of the Iterated Local Search is carried out through a Local Search strategy. The move of the Local Search consists of moving a single element from the allocation of one robot to another. We define two moves: -moveI (i, r 1 , r 2 ) : Note that moving a single element from one robot to another does not guarantee to produce an improvement in the objective function. Instead, to observe an improvement, it is frequently necessary to concatenate multiple moves. Thus, we accept all moves that do not decrease the objective. This strategy incurs the risk of making the Local Search rather ineffective and prone to cycling. To overcome this issue, we utilize heuristic information to guide the Local Search. Specifically, we estimate the benefit of all moves before attempting them so that only the most promising moves are considered.

Estimation of the move impact
The estimation of the move impact is carried out in three steps. First, we gather a set of local features of the scheduling, describing how each item and container relates to the rest of the current solution. Afterward, the gathered features are combined to describe a more general set of aggregated features, describing the status of the system at each point of the scheduling. Lastly, the features are evaluated to give an estimate of the quality of each available move. After having estimated the impact of the moves, we use a random approach to determine which move to be attempted, coupled with a short tabu list to prevent the method from cycling.

Computation of the local features
The first step of the evaluation computes the local features for the items and containers. All features are computed in relation to the current solution s and associated assignment A and scheduling T A . These features account for the contribution of the individual items and containers to the overall workload of the robots. In the following discussion, for the sake of brevity, we will only present the local features associated with the items.
The information we gather is summarized by a set of binary features for each robot-column pair (r , i), r ∈ R, i ∈ I with at least one item allocated, i.e., I r [i] ≥ 1. [i], let us consider the neighboring solution s, which corresponds to solution s with one more item i assigned to robot r . In this case, the addition of another item might have one of two effects. Either the added item i is discarded, and the decoded solution ofŝ corresponds with the decoded solution of s, or the item is not discarded, and the two decoded solutions differ. If the item is discarded, w r 6 [i] equals zero, otherwise it equals one. In the case of w r 7 [i], let us consider the neighboring solutionŝ, which corresponds to solution s with one less item i assigned to robot r . In this case, the removal of one item can have several possible outcomes. If an item i was discarded in the decoding of s, the decoded solutions ofŝ and s correspond and w r 7 [i] equals one. If no item i was discarded, the last item i associated to robot r in s is part of a mission (i, j, r ) p . Inŝ, slot j is either lost, or part of a different mission (i , j, r ) with i > i. If slot j is lost, w r 7 [i] equals zero. If slot j is part of a different mission (i , j, r ) p we further distinguish two cases: either t p (i, j,r ) < t p (i , j,r ) or t p (i, j,r ) = t p (i , j,r ) . In the first case w r 7 [i] equals zero, otherwise it equals one.
The local featuresw r 1−7 [i] of the containers are computed in a similar way.

Computation of the aggregated features
We now describe how the heuristic features discussed earlier can be aggregated into higher-level features, which describe the condition of the scheduling from a higher-level point of view. Those aggregated features are computed propagating the information obtained at the local level, to neighboring items, when it can be applied. They are evaluated iteratively, starting from the last items on the conveyor. For the purposes of the computation, we consider W r

Evaluation of the moves
After the features of the items and slots have been computed, the impact of moveI (i, r 1 , r 2 ) or moveC(i, r 1 , r 2 ) on the solution value is estimated. Let us define V r 1 − [i] and V r 2 + [i] as the heuristic value assigned to the removal of item i from robot r 1 and the addition of an item i to robot r 2, respectively. They are computed as follows: where l1, l2, g1, g2, g3, and g4 are positive constants. Formally, the heuristic evaluation of MoveI (i, r 1, r 2) is

Move selection scheme
To select the moves to be applied in the search, we systematically apply the move evaluation scheme presented in the previous section. The selection is performed iterating over the columns i ∈ I . For each iteration step, we compute the heuristic evaluation of all moves MoveI (i, r 1, r 2), for all r 1, r 2 ∈ R. A move is attempted if its evaluation V r 1,r 2 [i] exceeds a move selection threshold th. If none of the moves is attempted, we decrease the threshold by one and iterate to the next column. At the start of the search and whenever a move is attempted, the move selection threshold th is initialized to a random value in the interval [0, th max ]. If the attempted move was accepted, th is also increased by the current iteration number. This is done to encourage the method to explore moves uniformly.
To discourage cycling of the moves, we also implement a simple tabu mechanism in the search. When a move moveI (i, r 1, r 2) is successful, for the next tabu_length iterations of the search we mark as tabu any move where an item i is added in the scheduling of robot r 1, or where an item i is removed from the scheduling of robot r 2. On the other hand, if a move is rejected, we mark as tabu any move where an item i is added in the scheduling of robot r 2, or where an item i is removed from the scheduling of robot r 1. The moves that are marked as tabu are penalized by a factor of tabu_ penalt y in the heuristic evaluation of the move. The approach is summarized in algorithm 5.
Algorithm 5: The move selection process to select and apply the next moveI to the solution.

return;
The Local Search is stopped after N L S iterations without an improvement in the objective function.

Initial solution
The initial solution of the ILS is provided by a greedy constructive heuristic. Starting from an empty allocation, we construct the solution iteratively, building the solution of one robot at a time. For each robot, we generate the optimal single robot solution, with the remaining item and containers that are not part of the scheduling.
Practically, this is achieved starting from an empty allocation. For each robot, we temporarily allocate to it all the currently unallocated items and slots, compute the missions of the allocation, and remove all items and slots that are not assigned to any mission. The procedure is iterated until all robots have been processed.

Perturbation phase
The perturbation phase is used to escape the local optima reached at the end of the local improvement phase. To perturb the solution, we utilize a strategy to destroy a portion of the current solution, surrounding a critical item or container. Specifically, we select an item or a slot that is not assigned to any mission in the current solution. We remove from the allocation of the robots the selected element (column i ∈ I or container j ∈ C) and all items i ∈ I or slots j ∈ C such that |i − i| ≤ P I | j − j| ≤ P C , where P I and P C are parameters of the search.
We evaluate the perturbed solution and remove from the allocation all items and containers that did not fit in any mission. We then reconstruct the solution by applying the constructive heuristic presented in Sect. 5.5. To introduce additional diversification in the procedure, we reconstruct the allocation of the robots in a random order, instead of proceeding sequentially.
To explore uniformly the solution space, in our approach, we will alternate between applying the perturbation to items and to containers.

Extensions of the ILS
The metaheuristic approach we described so far is only able to handle the fixed speed version of the scheduling problem. In the next sections, we will discuss how to extend the metaheuristic method to handle the variable speed case and the item holding case.

Variable speed
We will use the same representation we proposed for the exact methods in Sect. 4.2.1. In addition, we define the following move as far as the speed is concerned.
In case the move reaches the minimum or maximum speed of the conveyor, we will limit the corresponding τ to that value. To avoid making moves with no effect on the speed profile, when this happens, we increase h by one.
The search for an effective speed profile for the problem is conducted using an iterative improvement strategy. At the start of the algorithm, we initialize the speed profile to a uniform speed across the planning horizon, with the speed set to the steady-state speed of the line. We then use a Local Search strategy to improve the initial speed profile.
To evaluate the quality of a speed profile, a high-quality assignment for that speed profile is required. Therefore, at each iteration we greedily generate an allocation for the line, and improve the assignment using the local improvement strategy presented in Sect. 5.4. We use the generated allocation to evaluate the quality of the current conveyor speed profile. All moves that do not deteriorate the solution are accepted.
We select the moves attempted in the search according to a greedy randomized heuristic. Let s be the incumbent speed profile. Additionally, let A be the assignment that has been generated for s. On A, we select a random lost item i and a random lost slot j. We note that if the current assignment does not have any lost item or slot, the current generated solution is already optimal, and the search may terminate.
Depending on the selection of i and j, two cases may occur, either a r i < α r j or a r i ≥ α r j , for a given reference robot r ∈ R. In the first case, we can hypothesize that the speed of the output conveyor may be increased, so as to move closer the availability intervals of i and j. In the second case, the speed of the output conveyor may be decreased to obtain the same effect. Therefore, in the first case we apply move inc (x, −μ, h), where x ∈ {1, . . . , k} : x j=1 τ j < a r i . In the second case, we apply move inc (x, μ, h), where x ∈ {1, . . . , k} : x j=1 τ j < α r j . The Local Search strategy is stopped after a total of N S1 + N S2 iterations. During the first N S1 iterations, we set μ = μ 1 and h = h 1 . Afterward, we set μ = μ 2 and h = h 2 . With μ 1 < μ 2 and h 1 > h 2 , doing so causes the initial moves to apply broad corrections to the speed profile and move the solution away from the steady-state solution, and later focus on finer adjustments of the profile. After N S1 + N S2 total iterations, we fix the speed profile and apply the ILS algorithm to obtain a final solution.

Item holding
The developed ILS can be easily adapted to allow for item holding. The adaptation does not require any structural change of the heuristic. The only change needed is to adjust the solution decoding procedure to account for the variable mission time.

Computational results
We performed a series of computational tests to evaluate the performance of the proposed approach and compare it to the more direct modeling approach used in Ferrari et al. (2015). The authors of Ferrari et al. (2015) provided us with data from three 4-manipulators plant configurations, representative of the characteristics of real-world lines where a dynamic scheduling approach could be applied. We report in Table 1 the parameters of the three considered lines are details. For each line, the table reports the column distance(d i ) and the maximum number of items in each column(n i ), the container distance(d j ) and the number of slots in each container(m j ), the velocity of    For each line configuration, we generated three instances with 15, 30, 45, and 60 columns of items, simulating the start-up conditions of the systems. Additionally, each of the generated instances was used as the basis for ten additional instances where some of the items are missing from the input conveyor and some slots are already filled at the start of the planning horizon.
The exact methods presented in this section were implemented with the AMPL and solved using CPLEX on a 2.00GHz Intel Xeon Processor L5335. The CPLEX solver uses 1 core, and the time limit has been set to 250 seconds. The metaheuristic developed have been written in C++ and ran on a 3.20Ghz AMD Ryzen 5 1600, in single-thread mode. The tuning parameters of the Iterated Local Search have been determined experimentally (Schettini 2017) and are reported in Table 9, in the "Appendix."

Exact methods
In Table 2, we report the summary of the results of the row generation approach in the basic version of the problem (base) and using the two presented extensions (holding, one speed, n speeds). For each line configuration and instance size, we report the number of instances solved to optimality (#S), the average number of lost items(obj), and the average solution time(t(s)). Note that the developed row generation algorithm does not provide a feasible solution when terminated early. Indeed, the developed cut generation algorithm solves a relaxed version of the problem and iteratively adds violated inequalities until a feasible solution is found. When the optimal solution to the relaxed problem is feasible, that solution is optimal by construction. Conversely, when the instance is not solved to optimality due to early termination, the solution of the relaxed problem is still not feasible; therefore, the algorithm only provides lower bound. Thus, for sake of fairness, we report its average objective only when all averaged instances are solved to optimality; when this is not the case, the value is reported as "-".
In the controllable conveyor case, two versions of the problem are considered, the case where the speed of the output conveyor can be controlled but can assume one value for the entire scheduling, and the case of a more finely controllable speed profile.
In practical terms, the first case (one speed) corresponds to having a speed profile composed of exactly one section, i.e., l = ∞. In the second case (n speeds), we consider sections of unitary length, i.e., l = 1. All versions of the problem were comparable in terms of computational time and instances solved.
Both extensions of the problem caused an improvement of the objective function. However, the benefits of the controllable conveyor strategy significantly outperform the holding strategy. The resolution of the speed profile did not affect the quality of the computed solutions, nor the difficulty of the problems, despite the increased freedom of the system.
The full experimental results are provided in the "Appendix" on Tables 6 and 7. Those tables, for each instance, report the best known lower bound of the optimal solution for each version of the problem, the computational time used to compute the solution, and the optimality of the solution ( * denotes solutions solved to optimality).

Comparison with Ferrari et al. (2015)
We compared our exact algorithm (ROW) with the model developed by Ferrari et al. (2015) (FERR) This comparison was carried out on the unmodified problem (base) and on the extended problem with controllable conveyor (one speed). In the latter case, to align with the model presented in Ferrari et al. (2015) we consider a constant speed velocity profile, meaning that the velocity profile is controllable, but only one speed can be selected for the entire planning horizon. Table 3 reports the summary of results. For each line configuration and instance size, we report the number of instances solved to optimality (#S), the average number of lost items(obj), and the average solution time(t(s)). The results are reported both for the basic problem and the controllable conveyors extension. The gap reported in the table is computed from the complement of the objective function, which maximizes the number of items moved into the containers. Table 5 reports the complete experimental results, for each instance, and for each solution method, the table reports the objective value and the solution time. We also report the optimality of the solution in the case of ROW and the optimality gap in the case of FERR.
We observe that ROW outperforms FERR in terms of computation time, being able to solve to optimality a larger portion of instances. However, we note that due to its nature, when the row generation does not reach an optimal solution, the solution provided at the end of the computation is not feasible. Conversely, the final solution provided by FERR is feasible, even if suboptimal. In the case of ROW, it is possible to repair the solutions obtained at the end of the computation utilizing the decoding algorithm presented in Sect. 5.1. However, on such occasions, it is more convenient to directly apply the ILS heuristic on the instance.
Overall, neither ROW nor FERR provide computational times that are compatible with a real-time application of the algorithm on any of the lines considered. However, their results are still useful for the purposes of understanding the problem and evaluating the developed heuristic.

Metaheuristic approach
In Table 4, we report the summary of the results of the developed heuristic approach on the 60 column instances. For each line configuration, we report the average number of items lost(obj), the best known lower bound to the number of lost items(LB), and the average solution time of the metaheuristic approach(t(s)). The problem has been solved in its basic version and considering the item holding and controllable conveyors extension. From the results, we observe that the computation time associated with the developed heuristic is limited to a couple of seconds, maintaining a high solution quality. Thus, the computational time achieved by the ILS is compatible with a real-time application of the algorithm to the considered lines. The full experimental results are provided in the "Appendix" in Table 8. The table, for each instance, reports the solution obtained by the Iterated Local Search, and the computational time utilized. We also report the best known lower bound of the corresponding instances.

Conclusions and future work
In this paper, we discussed the Pick and Place Packaging Problem for the optimization of a two-conveyor packaging line. We presented an integer linear programming model for the problem and an efficient exact algorithm based on row generation. We also developed an effective and time efficient heuristic based on Iterated Local Search. We presented two extensions of the problem. Specifically, we consider the case where one of the conveyors can be controlled, and the case where the robots may hold items without immediately mov-ing to their destination. Both the exact algorithm and the Iterated Local Search have been extended to handle these two extensions. The developed model and algorithms have been tested on a series of testing instances simulating realistic plant configurations. Our results show the effectiveness of the developed approaches and the utility of the optimization of robot scheduling. Some extensions of the problem can deserve further investigations. One is the multi-grip case when robots can pick up more than one item at a time and place them in multiple containers. Another case considers items of different types mixed in the system. In both cases, the extension of the mathematical model and of the heuristic can be particularly interesting.
Funding Open access funding provided by Politecnico di Milano within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.