Paper deep dive
Integrated packing, placement, scheduling, and routing of personalized production: a pharmaceutical Industry 4.0 use-case with a planar transport system
Viktor Emil Korladinov, Antonin Novak, Zdeněk Hanzálek, Erik Sonntag, František Štěpánek
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 98%
Last extracted: 4/27/2026, 4:29:53 AM
Summary
The paper proposes a multi-stage optimization framework for a planar Flexible Manufacturing System (FMS) applied to personalized pharmaceutical production. The framework decouples tactical decisions (long-term layout and resource packing) from operational decisions (daily scheduling and routing). At the tactical level, a Mixed-Integer Quadratic Programming (MIQP) model is used for drug packing to exploit co-occurrence patterns, followed by a placement problem to minimize expected travel distances. At the operational level, a Constraint Programming (CP) model treats movers as 'take-give' reservoir resources for scheduling, while an iterative conflict-resolution mechanism and DAG-based reasoning are used for conflict-free routing. The system is evaluated using real-world prescription data, demonstrating scalability and efficiency for personalized drug assembly.
Entities (8)
Relation Signals (5)
Planar Transport System → isatypeof → Flexible Manufacturing System
confidence 100% · The recent emergence of planar transport systems necessitates re-evaluation of Flexible Manufacturing Systems (FMS)
Mover → ismodeledas → Take-give Reservoir Resource
confidence 100% · modeling movers as reservoir resources... modeling them as take-give resources
Mixed-Integer Quadratic Programming → solves → Packing Problem
confidence 100% · A Mixed-Integer Quadratic Programming model is utilized for the packing problem
Constraint Programming → solves → Scheduling Problem
confidence 100% · This scheduling problem is formulated using Constraint Programming
Lékobot → uses → Planar Transport System
confidence 90% · the “Lékobot” platform has been proposed... To achieve industrial throughput, the system employs planar transport
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:The recent emergence of planar transport systems necessitates re-evaluation of Flexible Manufacturing Systems (FMS) to address the simultaneous scheduling of internal logistics and production operations. By operating on a tile-based planar grid, these systems allow independent movers full two-dimensional freedom, mitigating inefficiencies inherent to traditional sequential lines. This paper applies a planar FMS framework to a real-world use case in the pharmaceutical industry: the automated production of personalized drugs. Implementing this system requires solving optimization problems at both tactical and operational levels. The tactical level involves decisions regarding production line layout and the positioning of drug dispensers. A Mixed-Integer Quadratic Programming model is utilized for the packing problem to exploit drug co-occurrence patterns found in historical patient data. Subsequently, we solve the placement problem - a bi-level problem combining an assignment problem with Shortest Hamiltonian paths with neighborhoods - to arrange dispensers in a layout minimizing expected travel distances. The operational level is encountered daily, scheduling individual movers to process new orders as quickly as possible. This scheduling problem is formulated using Constraint Programming, modeling movers as reservoir resources to ensure order completeness, complemented by a routing phase using an iterative conflict-resolution mechanism and DAG-based reasoning to convert schedules into conflict-free paths. Evaluation using real-world prescription data for 40 drugs shows the framework scales efficiently across several layout topologies for up to 500 orders, with schedules that are highly effective and computationally tractable for daily operations.
Tags
Links
- Source: https://arxiv.org/abs/2604.21029v1
- Canonical: https://arxiv.org/abs/2604.21029v1
Trouble viewing inline? Open PDF directly →
Full Text
100,117 characters extracted from source content.
Expand or collapse full text
Integrated packing, placement, scheduling, and routing of personalized production: a pharmaceutical Industry 4.0 use-case with a planar transport system Viktor Emil Korladinov a,b,∗ , Antonín Novák a , Zdeněk Hanzálek a , Erik Sonntag c , František Štěpánek c a Czech Institute of Informatics, Robotics and Cybernetics, Czech Technical University in Prague, Jugoslávských partyzánů 1580/3, 160 00, Prague 6, CZ b Faculty of Electrical Engineering, Czech Technical University in Prague, Karlovo Náměstí 13, 120 00, Prague 2, CZ c Department of Chemical Engineering, University of Chemistry and Technology, Technická 3, 166 28, Prague, CZ Abstract The recent emergence of planar transport systems necessitates re-evaluation of Flexible Manufacturing Systems (FMS) to address the simultaneous scheduling of internal logistics and production operations. By operating on a tile-based planar grid, these systems allow independent movers full two-dimensional freedom, mitigating inefficiencies inherent to traditional sequential lines. This paper applies a planar FMS framework to a real-world use case in the pharmaceutical industry: the automated production of personalized drugs. Implementing this system requires solving optimization problems at both tactical and operational levels. The tactical level involves decisions regarding production line layout and the positioning of drug dispensers. A Mixed-Integer Quadratic Programming model is utilized for the packing problem to exploit drug co- occurrence patterns found in historical patient data. Subsequently, we solve the placement problem — a bi-level problem combining an assignment problem with Shortest Hamiltonian paths with neighborhoods — to arrange dispensers in a layout minimizing expected travel distances. The operational level is encountered daily, scheduling individual movers to process new orders as quickly as possible. This scheduling problem is formulated using Constraint Programming, modeling movers as reservoir resources to ensure order completeness, complemented by a routing phase using an iterative conflict- resolution mechanism and DAG-based reasoning to convert schedules into conflict-free paths. Evaluation using real-world prescription data for 40 drugs shows the framework scales efficiently across several layout topologies for up to 500 orders, with schedules that are highly effective and computationally tractable for daily operations. 1. Introduction The manufacturing sector is currently undergoing a structural transformation driven by the demand for mass personalization. Unlike traditional mass production, which relies on rigid, linear workflows to achieve efficiency, modern production environments must be able to handle diverse manufacturing procedures with ∗ Corresponding author Email addresses: korlavik@fel.cvut.cz (Viktor Emil Korladinov), antonin.novak@cvut.cz (Antonín Novák), zdenek.hanzalek@cvut.cz (Zdeněk Hanzálek), Erik.Sonntag@vscht.cz (Erik Sonntag), Frantisek.Stepanek@vscht.cz (František Štěpánek) 1 arXiv:2604.21029v1 [math.OC] 22 Apr 2026 minimal setup and transportation times. Flexible Manufacturing Systems (FMS) have already been heavily studied as a solution, e.g., see survey paper by Yadav & Jayswal (2018). However, the recent emergence of planar transport systems (e.g., Beckhoff XPlanar 1 ) creates new opportunities and necessitates a re-evaluation of FMS as it represents a significant technological leap. While two-dimensional transport is not new to FMS, earlier vehicle-based systems usually required robotic arms’ assistance to interact with stations. Modern tile-based systems achieve sufficient positioning precision for movers to interact directly with dispensers, eliminating such overhead along with the inefficiencies inherent to sequential lines. While this hardware offers unprecedented flexibility, it introduces a massive combinatorial design space exploration that we tackle in this paper. Unlike traditional lines, where the sequence of operations is dictated by the layout, a planar grid allows for independent routing and station sequencing for each order. This creates a complex multi-stage optimization challenge that must integrate tactical decisions (such as resource packing and spatial placement) with operational decisions (including scheduling and conflict-free routing). As shown by Maenhout & Vanhoucke (2013), decoupling these tactical decisions from daily operations yields good results and robust performance, provided the tactical plans anticipate operational constraints. However, existing frameworks often treat these problems in isolation. Standard Multi-Agent Path Find- ing (MAPF) approaches (reviewed in depth by Stern (2019)) typically prioritize collision avoidance for fixed layouts, ignoring the broader impact of station placement on scheduling decisions. Engineering intuition is insufficient for such coupled systems, as minor adjustments in station placement can drastically alter throughput. Consequently, there is a critical need for formal optimization frameworks that can mathemati- cally model these interactions and determine cost-effective configurations prior to physical implementation. A simulation tool to help visualize the framework’s work and output is available on YouTube 2 . 1.1. Motivation: Personalized Pharmaceutical Production We apply this planar FMS framework to a real-world use case in the pharmaceutical industry: the au- tomated production of personalized medication. The efficient management of polypharmacy is a critical healthcare challenge, with an increasing number of patients prescribed five or more concurrent medica- tions (Wastesson et al., 2018). Although the pharmaceutical landscape is surprisingly narrow (39 medica- tions account for approximately 50% of all prescriptions (Quinn & Shah, 2017)), the vast number of potential drug combinations and dosage strengths creates immense combinatorial variability. This variability renders the traditional mass-production model of Fixed-Dose Combinations economically unfeasible. As noted by Hao et al. (2015) and Sutherland et al. (2015), regimens containing more than three active ingredients become effectively unique to the patient, leading to excessive Stock Keeping Unit proliferation if manufactured via traditional means. Current solutions are polarized between rigid industrial platforms, which lack flexibility (Vanhoorne & Vervaet, 2020), and manual pharmacy compounding, which lacks scalability (Carvalho & Almeida, 2022). To bridge this gap, the “Lékobot” platform has been proposed as a modular approach in which individ- ualized doses are assembled from prefabricated subunits (mini-tablets or pellets) into capsules (illustrated in Figure 1c, motivated by Sonntag et al. (2023)). This moves dose adjustment from the formulation level to an automated assembly process. To achieve industrial throughput, the system employs planar transport 1 Please see XPlanar system showcase (≈ 1min long) https://w.youtube.com/watch?v=yIiZ_DzV0c8&t=60s 2 Simulation tool showcase (≈ 2min long) https://w.youtube.com/watch?v=ukVYjog60Js 2 (Figure 1d) to route patient-specific batches through a series of dispensers (shown in Figure 1a). This ap- plication provides a rich instance of the planar FMS problem, requiring the rigorous optimization of layout, resource packing, routing, and scheduling to balance the flexibility of compounding with the efficiency of industrial manufacturing. 1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 A B C D E F (a) A tile, a mover carrying a 100-capsule cartridge, and an overhead dispenser. (b) Robot arm interface designed for the loading and unloading of capsule cartridges. (c) Capsules filled with mini-tablets. 1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 A B C D E F (d) Planar system with 25 tiles and 4 movers. Figure 1: Components of personalized drug production. 1.2. Contributions and Paper Outline The main contributions of this paper are as follows: • We formalize the automated assembly of personalized medication on a planar transport system. • We introduce a comprehensive, multi-stage optimization framework that decouples tactical, long-term decisions (resource packing and spatial placement) from operational, daily execution (scheduling and routing). • At the tactical level, we formulate a Mixed-Integer Quadratic Programming (MIQP) model to exploit drug co-occurrence patterns, and we apply a Genetic Algorithm (GA) with the inner Shortest Hamil- tonian path with neighborhood problem to optimize the physical placement of dispensers to minimize expected travel distances. • At the operational level, we develop a Constraint Programming (CP) model that treats movers as take-give reservoir resources, coupled with a Directed Acyclic Graph (DAG) based routing logic to obtain collision-free paths. 3 • We prove the complexity of the placement problem and introduce a lower bound on the makespan that combines parallel machine scheduling with the Noon-Bean transform. • We evaluate the proposed approach using synthetic data derived from real-world prescription patterns, demonstrating that our framework consistently produces near-optimal makespan while effectively scal- ing through a batch-merging heuristic. The remainder of this paper is organized as follows. Section 2 formally defines the problem at both the tactical and operational levels. Section 3 reviews related work in flexible manufacturing and multi-agent pathfinding. The tactical level is addressed in Section 4 (Packing) and Section 5 (Placement). The operational level is covered in Section 6 (Scheduling) and Section 7 (Routing). Section 8 presents the experimental setup and discusses the results. Finally, Section 9 concludes the paper. 2. Problem Statement As described above, we are faced with several optimization problems at two different levels. The tactical level is where long-term decisions, such as the layout of the production line and the number and positions of dispensers and interfaces, are made. For example, the production line will be constructed in a certain way and will remain in its current configuration for many months. The second level—the operational, is encountered on a daily basis when a new set of orders is received, with the goal of scheduling individual movers to produce it as soon as possible. A diagram illustrating the relationships between the two decision levels is shown in Figure 2. Tactical level Operational level dispensers n disp maximum dispensers d max historical orders P history tiles n tiles maximum movers m max Packing 01234567 tile utilization layout H, |H| = n tiles + n inter interfaces I, |I| = n inter Placement new customer orders P available movers M, |M| = n movers SchedulingRouting packing D placement φ :D∪I →H schedules m 1 m 2 m 3 Figure 2: Diagram of the proposed solution and mutual interactions of its components. 4 2.1. Tactical Level The goal of the tactical level is to make long-term decisions related to the hardware construction of an efficient production line. For this, we have access to a set of historical orders P history that were processed in the past, which we assume to be a representative sample for making tactical decisions. Each observation in the dataset represents an individual participant, identified by a unique sequence number. Each participant has a list of prescribed pharmacological agents and their corresponding dosage frequencies. Namely, at the tactical level, we determine how many dispensers we should allocate to each drug and assign these dispenser locations within the given layout of the transportation system on the factory shop floor. The layout is a subset of 2D integer lattice H ⊆N×N. Each element (i,j) ∈ H corresponds to an available tile in the system at (x,y) coordinates. For example, a square layout with total of 64 tiles as shown in Figure 5a, is represented as H square 64 = 1,..., 8 2 and the ring layout with 64 tiles as shown in Figure 5b is given as H ring 64 = 1,..., 17 2 \2,..., 16 2 . The tactical-level problem utilizes a set of n disp dispensers to be placed over n tiles dispensing tiles to provide the supply of a set of drugs G by using the maximum number m max of available movers. The mover starts each order on an arbitrary interface tile (to receive an empty cartridge) and ends on an arbitrary interface tile (to release a filled cartridge; the two interfaces can, but do not have to be the same). There are n inter =|I| interface tiles. Each dispenser contains one type of drug (to avoid contamination), but any drug k ∈G can be allocated to multiple dispensers. Further, we assume that each tile can accommodate at most d max dispensers due to space constraints (typically 4) or one interface location i∈I. The upper green area in Figure 2 represents the tactical level of decisions. In this work, we make a deliberate choice to split the tactical level decisions into two problems. First, in Section 4, we derive a Packing problem that (i) decides how many dispensers should be allocated for each drug and (i) packs specific dispensers into unplaced dispensing tiles to leverage typical drug co-occurrence in individual orders. At this stage, a dispensing tile is unplaced—it only represents a logical grouping of drugs that will eventually share a dispensing tile together, though the physical location of that dispensing tile on the grid remains undecided. Formally, the Packing problem associates drugs with these unplaced dispensing tiles to create a set D =(g,h)| g ∈G,h∈1,...,n tiles denoting that drug g is assigned to a specific (yet unplaced) dispensing tile h. The subsequent Placement step translates these logical groupings into a physical configuration. Having defined the composition of each tile, we must now determine their specific coordinates within the production line. As detailed in Section 5, this involves mapping both the pre-packed dispensing tilesD and the interfaces I =(interface,h)| h∈1,...,n inter onto the physical layout H. We seek an assignment φ :D∪I →H where each available coordinate (x,y)∈H is occupied by either a single interface i ∈ I or a dispensing tile. In the latter case, the tile at (x,y) carries the specific set of dispensersg | (g,h ′ )∈D previously grouped into unplaced dispensing tile h ′ during the packing stage. By optimizing this spatial arrangement, we aim to minimize the makespan of operational schedules for future orders P, assuming they follow the same distribution as the historical data P history . 2.2. Operational Level The operational level provides decision-making support to run daily operations, ensuring that orders planned for each day are produced. Hence, the goal of the operational level is to determine how to operate the number of available movers n movers on the given layout of dispensers φ so that we produce the set 5 of current orders P. Each order p ∈ P is given by a set of all drugs S p ⊆ G. The required dispensing quantity of the drug g in order p is an integer ∆ g p ∈N. Hence, the order p is represented as a set of tuples p =(g, ∆ g p )|∀g ∈S p . To complete an order p, the mover starts and stays for η I ∈N time units at one of the interface locations (to mount an empty cartridge with capsules), visits all dispensers with required drugs g ∈S p in an arbitrary order (i.e., in open-shop manner), spends the requested dispensing time ∆ g p there, and finishes again in an interface location, spending again η I time units. Thus, at any moment, the mover processes only one order at a time, and each interface or dispenser tile serves a single mover. The dispensing process can be preempted at any time and later resumed. The paths of the movers are Manhattan-style (i.e., follow ℓ 1 distance) and are conflict-free. The goal is to minimize the makespan, that is, the latest completion time of any mover. See an illustration in Figure 2 where the bottom green area displays all required inputs for making decisions at the operational level. Similarly to the tactical problem, we deliberately decide to split the operational problem into two sub- problems here. First, in Section 6, we deduce the Scheduling problem that assigns orders p ∈ P to the available amount of moves M and decides on the sequencing of the individual dispensing operations. The setup times between dispensing operations, reflecting the distance between dispensers, are enforced, but we choose to relax on mover conflicts during their transports. To handle these conflicts, we introduce in Section 7 the Routing problem that finally produces conflict-free paths for the movers by slightly extending the optimistic schedule generated from the Scheduling problem solver. 3. Related Work Flexible Manufacturing Systems (FMS) are well-established in the literature, typically characterized as medium-volume automated systems (Yadav & Jayswal, 2018). However, the emergence of high-precision planar transport systems, such as Beckhoff XPlanar or B&R ACOPOS 6D, combined with the Industry 4.0 shift toward mass personalization (Wang et al., 2017), has increased system flexibility to a point that requires updated planning strategies. Motivated by these advancements, this work revisits flexible systems, distinguishing between the tactical and operational decision levels required for effective management. Maenhout & Vanhoucke (2013) demonstrate that decoupling tactical decisions from operational deci- sions is essential for robust performance, provided that tactical plans anticipate operational constraints. In the context of manufacturing, this allows us to optimize the layout and dispenser configuration (tactical) separately from the daily routing of movers (operational), similar to the robust recovery policies proposed by Akbarzadeh & Maenhout (2024) for healthcare resources. The operational level handles day-to-day order fulfillment. The need to assign operations to one of several available dispensers mirrors the routing flexibility of the Flexible Job Shop Problem, as reviewed by Dauzère-Pérès et al. (2024). The drug dispensing process allows for free operation ordering, a core feature of the Open Shop Problem (Abreu & Nagano, 2024). The coordination of movers requires modeling them as renewable resources. Hartmann & Briskorn (2022) classify such complex constraints in their survey of Resource-constrained Project Scheduling Problem extensions. What distinguishes this system from classical formulations is its reliance on internal logistics rather than machine reconfiguration. Traditional manufacturing lines are hindered by “long line” inefficiency; products must traverse every station regardless of necessity, as seen in customization systems like the “Lip Picker” (Go, 6 2021). Planar systems utilize dynamic routing to visit only required dispensers. This modularity allows the system to be reconfigured on the fly to manufacture different drugs. However, this transport-driven flexibility introduces a complex scheduling constraint: the coordination of the transport units themselves. Laborie (2003) defines a reservoir resource, which functions like a fuel tank: it has a capacity and is only assignable when that capacity is positive. A specific subclass of this is the take-give resource proposed by Hanzálek & Šucha (2017). In their model, capacity is binary (0 or 1); the resource is "taken" at the start of a task chain and released ("given") only upon completion. Hartmann & Briskorn (2022) gives structure to the different types of resources and classifies take-give resources as partially renewable cumulative resources. In this work, we model the movers as take-give resources, which ensures they remain dedicated to a single order until the entire dynamically routed sequence is complete. Despite the complexity of the scheduling logic, the physical execution on the planar transport grid allows for some routing flexibility: with our assumed XPlanar tile configuration, two movers can share a single tile, weakening the strict vertex-collision model assumed by classical Multi-Agent Path Finding (MAPF) (Stern, 2019). Even purpose-built approaches for levitating planar systems (Nilsson & Ternerot, 2022; Rehme, 2025) address only point-to-point navigation, and while Janning et al. (2025) successfully combine CBS with MPC to guarantee collision-free execution, none of these methods account for the multi-stop paths inherent to complete order fulfillment. Decoupled formulations such as Zahradka et al. (2022) separate task assignment from routing in the Multi-Agent Multi-Item Pickup and Delivery setting, but applying exact MAPF solvers within an iterative scheduler remains computationally prohibitive: as Fioravantes et al. (2024) demonstrate, MAPF is theoretically intractable and NP-hard even in basic planar graphs with fixed makespan. We therefore approximate routing logistics via sequence-dependent setup times and the take-give resource mechanism, augmented by an iterative conflict-resolution step, omitting MAPF frameworks from our solution. 4. Packing Problem The goal of the Packing problem is to decide how many dispensers will be allocated to each drug, and also which dispensers will be placed together on the same tile to take advantage of typical drug co-occurrences and to use the tile efficiently for less frequent drugs, as shown in Figure 2. The input consists of a list of historical orders P history (same structure as P), the total number of tiles n tiles ∈N, the total number of dispensers available n disp ∈N, the maximum number of dispensers on a single tile d max ∈N and maximum number of movers m max ∈N. The maximum number of dispensers to be placed on a tile is a parameter that follows from physical constraints imposed by the dimensions of dispensers and tile size. We remind that the goal of the packing phase is not to assign tiles into the layout, but only to decide on the (i) multiplicity of each drug and (i) co-occurrence of drugs on the tiles. Therefore, the problem solved by the packing phase is formulated at the capacity level, whereP history , the list of past orders, is used to estimate the future demand for each type of drug. The problem is formulated as the following MIQP model. We are given the following parameters: n tiles ∈N number of tiles available, n disp ∈N number of dispensers, u g ∈R + 0 estimated total cumulative dispensing time for drug g ∈G. This number essentially represents, e.g., the total time in seconds that will be spent by dispensing drug g ∈ G over all orders in P history . 7 As for the decision variables, we use y k,g ∈ 0, 1 indicator whether drug g is assigned on unplaced tile k (i.e., which is not placed at any particular position at the moment), z g ∈N 0 is amount of dispensers for drug g, π g ∈R + 0 is the estimated load for one dispenser for drug g, i.e., π g = u g /z g , assuming they will be utilized uniformly. Indeed, at the tactical stage of decision-making, we do not have access to the information about the scheduling of the specific movers (which is to be revealed at the operational level), hence we assume a uniform distribution of the load. Furthermore, μ k ∈R + 0 is the total load of all drugs on tile k and μ max = max k μ k is the maximum total load over all tiles. Since we have n inter to be placed on separate tiles without any dispensers, they are not considered in the model, thus the total amount of tiles is decreased by n inter . The full model is stated as follows: minμ max (4.1) subject to μ max ≥ μ k ∀k ∈1,...,n tiles (4.2) π g · z g = u g ∀g ∈G(4.3) n tiles X k=1 y k,g = z g ∀g ∈G(4.4) 1≤ X g∈G y k,g ≤ d max ∀k ∈1,...,n tiles (4.5) X g∈G π g · y k,g = μ k ∀k ∈1,...,n tiles (4.6) n tiles X k=1 X g∈G y k,g ≤ n disp (4.7) where y k,g ∈0, 1 ∀k ∈1,...,n tiles ,∀g ∈G(4.8) z g ∈0,...,m max ,π g ∈R + 0 ∀g ∈G(4.9) μ k ∈R + 0 ∀k ∈1,...,n tiles (4.10) μ max ≥ 0.(4.11) The output of the model is a tile-drug matrix, as shown, for example, in Figure 3a. Drugs are sorted along the horizontal axis in a non-increasing order of their expected total utilization u g . The available individual tiles are displayed along the vertical axis without any particular ordering, as they have not yet been placed in the layout. Each dot in the matrix corresponds to a value y k,g = 1, meaning that drug g ∈ G will be placed onto tile k. For example, we can see that the most requested drug will be allocated to 4 tiles. On the other hand, drugs along the same row share a common tile. The model objective (4.1) aims to minimize the maximum expected tile load, as shown in the solution displayed in Figure 4 and also by the shades of blue in Figure 3. Each bar in Figure 4a corresponds to the expected tile utilization μ k , which is calculated using equation (4.6), which is expressed as a sum of contributions of assigned dispensers on the tile k. The contribution of one dispenser π g for drug g ∈ G is calculated using the total demand u g for drug g, divided by the number of dispensers z g for drug g. This 8 010203040 50 40 30 20 10 0 drug g [-] tile k [-] (a) Tile-drug matrix of drugs using 65 dispensers. 010203040 50 40 30 20 10 0 drug g [-] 0 200 400 600 expected tile utilization μ k [-] (b) Tile-drug matrix of drugs using 100 dispensers. Figure 3: Tile-drug matrix represents values of variables y k,g for different amounts of dispensers n disp with n tiles = 50 tiles. is implemented by a bilinear constraint (4.3). Hence, we assume that if z g dispensers for drug g will be used, then the overall demand of g would be equally split among them, leading to their individual utilization u g /z g . This constraint expresses our assumption that, for high-quality placements, the operational-level schedule leads to balanced utilization of all tiles with drug g. The rest of the constraints ensure that we can fit at most d max dispensers onto a tile by constraint (4.5) and that the total amount of dispensers is not exceeded (4.7). Furthermore, we included the constraint (4.9) for practical reasons to limit the number of dispensers of drug g. After initial experiments, we observed that when n disp is large enough, then the model tends to generate many dispensers of the same drug to improve utilization balance, sometimes in such high numbers that they would be even larger than the movers available. 01020304050 200 400 600 tile [-] utilization μ k [-] (a) Expected tile utilization after packing using n disp = 65 dis- pensers, resulting in μ max = 732. 01020304050 200 400 600 tile [-] 1st dispenser 2nd dispenser 3rd dispenser 4th dispenser (b) Expected tile utilization after packing using n disp = 100 dis- pensers, resulting μ max = 661.5. Figure 4: Expected tile utilization with individual dispensers’ contributions. With the increasing number of dispensers n disp available, we can allow more balanced utilization of the tiles. For example, a solution to the packing problem with n disp = 100 is shown in Figure 3b with expected tile utilization plotted in Figure 4b—we see that maximum utilization is decreased compared to Figure 4a, as well as the tiles have a more balanced workload. As the number of tiles n tiles = 50 remains the same in both cases, more dispensers are being packed onto a single tile, as visible in Figure 3b. In such cases, it would be advantageous to assign dispensers to tiles 9 together, not only to minimize the tile workload but also to account for mutual drug co-occurrence, thereby reducing mover movements when completing an order. We have formulated the maximization of pairwise correlations as another optimization step following the minimization of the maximal utilization. Please see details of this model in the Appendix B. 5. Placement problem We assume that the input to the Placement problem consists of the union of the set of (unplaced) tiles D with packed dispensers together with the set of interface tiles I (used for extraction of cartridges from movers with completed orders and intersection of an empty cartridge for a new order), the given layout H, |H| = n tiles + n inter and a set of orders P history that reflect an empirical data, e.g., based on past orders. In this work, for simplicity, we assume four basic types of layouts: line (e.g., H line 25 = (1, 1),..., (1, 25)), doubleline, ring, and square, as shown in Figure 5. We refer to specific layout sizes with the string|H|∼|I|. For example, the square 8× 8∼ 4 layout refers to a square with 8 rows, 8 columns, and 4 interfaces. (a) Square 8× 8∼ 4 layout. (b) Ring 64∼ 4 layout. (c) Line 25∼ 1 layout. (d) Doubleline 2× 25∼ 2 layout. Figure 5: Examples of considered layout topologies. Highlighted tiles display an example of a possible placement of interface locations I, which is together with packed tiles D, subject to optimization. The goal of Placement is to find a mapping φ : D ∪ I → H between unplaced tiles with packed dispensers D (i.e., output of the Packing problem) and interfaces I to the layout H, such that the mean number of movers’ steps in Manhattan distance needed to fulfill a set of historical ordersP history is minimized: min φ:D∪I→H placement(φ) =min φ:D∪I→H 1 |P history | X p∈P history min σ∈I×Π(p)×I |σ|−1 X j=1 ∥φ(σ j )− φ(σ j+1 )∥ 1 ,(5.1) where I is the set of all interfaces and Π(p) is the set of all permutations of dispensers relevant to order p defined as follows. The set of dispensers providing drug g is denoted as D g =(g,k)|∀(g,k)∈D, i.e., the set of alternative dispensers of drug g. Then, Π(p) = Perm S (g,∆ g p )∈p D g be the set of all permutations of D g for all drugs g requested in the customer order p. Finally, ∀σ ∈ Π(p) : σj denotes j-th element of permutation σ of dispensers relevant to p. The Placement problem, as defined by (5.1), consists of two levels: (i) an assignment problem that assigns packed dispensers and interfacesD∪I to the specific tiles placed in the given layoutH, and (i) sum of the inner problems, where each resembles the Shortest Hamiltonian path problem with neighborhoods (SHPPN) (Arkin & Hassin, 1994; Deckerová et al., 2024). Each customer order p ∈ P history represents a single path in a graph defined by assignment φ, starting and finishing at some interface location i ∈ I. The path must contain the required dispensers, each of which may exist in multiple alternatives, forming 10 a so-called neighborhood. Due to the considered layouts and the permitted movements, the costs of the underlying graph resemble the ℓ 1 metric on an integer lattice, i.e., Manhattan distance. See an example in Figure 6a, where the solution φ of the assignment part of the Placement problem is displayed. To calculate the quality of placement φ, we have a set of three historical orders P history given in the drugs column of Figure 6b. With these, the optimal solution of the corresponding three inner problems in (5.1) is given in permutation σ column in Figure 6b. These correspond to the optimal sequence of dispensers together with a starting and finishing interface for each order in P history . Finally, the objective value of placement φ is calculated as the arithmetic average of the mover steps in each permutation σ. Hence, with the increasing size |P history | → ∞, the objective function approaches the expected number of steps needed to fulfill an order. The objective function of this bi-level problem is designed to find a placement φ such that dispensers with drugs requested in typical individual orders are placed near each other, thereby reducing the average travel time of movers. Later in the experiments in Section 8.3, we will show that the choice of the objective function (5.1) is justified by its correlation with the actual schedule makespan observed at the operational level. 1234 1 2 3 4 x coordinate y coordinate INTERFACE OMEPRAZOLE LEVOTHYROXINE LOVASTATIN VALSARTAN METFORMIN GLIPIZIDE PRAVASTATIN LISINOPRIL SIMVASTATIN METOPROLOL CLOPIDOGREL LOSARTAN AMLODIPINE ATORVASTATIN WARFARIN ATENOLOL FUROSEMIDE HYDROCHLOROTHIAZIDE (a) Example placement φ in square layout 4× 4∼ 2. order pdrugspermutation σsteps [-] #1atorvastatin, hydrochlorothiazide(2, 1)−→ (3, 1)−→ (3, 2)−→ (3, 3)3 #2omeprazole(3, 3)−→ (2, 3)−→ (1, 3)−→ (1, 4)−→ (2, 4)−→ (3, 4)−→ (3, 3)6 #3lisinopril, simvastatin(3, 3)−→ (2, 3)−→ (2, 2)−→ (2, 1)3 placement(φ) = (3 + 6 + 3)/3 = 4 (b) Example P history data and their optimal permutations σ. Figure 6: Example of the considered placement in the example square layout 4× 4∼ 2. 5.1. Complexity of the Placement problem In this section, we hint at two sources of its complexity. The first one comes from the assignment part of the problem, agnostic to the SHPPN arising in the inner problem. We show that the decision version of Placement is NP-hard even for the line layout. Proposition 1 (Complexity of Placement). Deciding whether Placement has a solution with the objective value less than or equal to q is NP-hard. The reduction will be shown from the Exact 3-Cover problem: Definition 1 (X3C). Exact 3-Cover: 11 • Given: Set U =1, 2,..., 3q, family of sets F =S 1 ,S 2 ,...,S k , |S i | = 3, S i ⊆ U. • Goal: Does exist an exact cover of U, i.e., exists F ′ ⊆ F, |F ′ | = q such that S i ,S j ∈ F ′ , i ̸= j: S i ∩ S j =∅ and S S i ∈F ′ S i = U. X3C is known to be NP-complete, even if each i ∈ U is contained in at most three sets S i (Garey & Johnson, 1979). Proof. We construct the following instance of Placement problem. Consider a line layout with k + 2 tiles containing 2 interfaces. The set of all drugs corresponds to the set U. The packing of dispensers is given by each S i ∈ F set, i.e., we have 3k dispensers packed into a total of k (unplaced) tiles. Finally, we have a single order in P history that requires all drugs U. We claim that X3C has a solution if and only if there is a placement φ with the cost equal to q mover steps. Suppose we have solution F ′ =S ′ 1 ,...,S ′ q to X3C. We construct a solution to the placement problem in the following way. Place the first interface to tile (1, 1) and the second interface to tile (1,q + 2). The q tiles that pack a set of dispensers defined by the solution of X3C S ′ 1 ,...,S ′ q are placed between the interfaces in an arbitrary order. The remaining k−q tiles are placed in arbitrary order on tiles (2,q+2), (1,q+3),..., (1,k+2). See visualization in Figure 7. S ′ 1 ,...,S ′ q Figure 7: Placement solution corresponding to a feasible instance of X3C. It can be seen that such a solution has cost q and satisfies the customer order requirement, since a mover fulfilling the order can traverse from the left to the right interface, dispensing all 3 drugs on each visited tile. Further, note that no better solution exists as the order requires|U| = 3q drugs, and each tile is packed with exactly 3 dispensers. Therefore, the mover has to visit at least q tiles. Next, consider a feasible solution to the placement problem with a cost equal to q, shown in Figure 8. Let a and b, a < b be indices of tiles where interface locations are placed. Further, let S ′ a+1 ,...,S ′ b−1 be the sets defined by drugs that are dispensed at tiles (1,a + 1),..., (1,b− 1). Since the order requires 3q = |U| drugs and each tile (except an interface) contains 3 dispensers, then no such solution can start and finish at the same interface location. Therefore, in the optimal placement, there are b− a− 1 = q tiles between the interfaces. Since any feasible solution of placement with cost q for a single order must ensure that all its 3q drugs can be dispensed in q steps of the mover, it follows ∀i,j ∈a + 1,...,b− 1,i < j : S ′ i ∩S ′ j =∅, hence F ′ =S ′ a+1 ,...,S ′ b−1 is an exact cover of U. S ′ a+1 ,...,S ′ b−1 Figure 8: A solution of placement problem where q steps of a mover fulfills the order drug requirements U. Note that when we allow the packing of at most 2 dispensers onto the same tile, the above reduction would no longer hold, as the Exact 2-Cover (covering by sets with exactly two elements) is polynomially solvable—can be reduced to 2-SAT, which is polynomially solvable (Garey & Johnson, 1979). Thus, this leaves the possibility of having lower complexity on simple topologies such as line or ring if the maximum number of dispensers on a tile is limited to 2. However, as indicated above, the problem of placing dispensers and interface locations is not the only source of hardness. Without providing a formal proof, we suspect the problem may still be hard due to the inner SHPPN problem. Indeed, it is known that the standard 12 Traveling Salesman Problem in the plane equipped with ℓ 1 metric still remains strongly NP-hard, see (Itai et al., 1982) or [ND23] in (Garey & Johnson, 1979). 5.2. Genetic Algorithm with Inner SHPPN Problem To solve the Placement problem represented by equation (5.1), we proposed the following algorithm. Due to the complexity of the inner problem, we proposed a meta-heuristic algorithm that assigns packed tiles D to the layout H, where the placement quality of φ is approximated as follows. For each customer order p in P history , we solve the inner optimization problem in the form of [SHPPN]min σ∈I×Π(p)×I |σ|−1 X j=1 ∥φ(σ j )− φ(σ j+1 )∥ 1 ,(5.2) which involves finding the shortest path from an interface location over a set of selected vertices representing relevant dispenser locations for order p, finishing in some interface location as well. The underlying graph resembles a lattice equipped with ℓ 1 metric. Hence, the problem effectively reduces to finding a sequence σ of dispensers that the mover visits, as the path between locations φ(σ j ) and φ(σ j+1 ) always exists with length exactly ∥φ(σ j )− φ(σ j+1 )∥ 1 . For example, an optimal solution to these inner problems is shown in Figure 6b for the placement from Figure 6a, for three historical customer orders. If one would not solve SHPPN problem exactly, one of the natural ways of solving (5.2) would be to find σ greedily—start at some interface location i ∈ I and consecutively proceed from φ(σ j ) to the following location φ(σ j+1 ) with the lowest distance ∥φ(σ j )− φ(σ j+1 )∥ 1 and so on. However, this obviously may lead to suboptimal results, and at the operational level, it may not always be possible, for example, because the nearest required dispenser may be occupied by another mover. Therefore, at the operational level, we would sometimes have to assign a more distant dispenser, i.e., breaking the greedy ordering. Hence, instead of a fully greedy method, we construct multiple σ, where each sequence of dispensers is also constructed in a greedy way, but the next dispenser in the sequence is not taken greedily, but is sampled inversely proportionally with respect to their distance from the current location φ(σ j ). Thus, closer alternative dispensers are more likely to be included in the calculation of (5.2) (reflecting what an efficient scheduler at the operational level would do), while not ruling out the occasional inclusion of more distant ones. Consider, for example, the optimal permutation for order #1 in Figure 6b. An alternative permutation σ ′ would be (2, 1) −→ (2, 2) −→ (2, 3) −→ (3, 3) −→ (4, 3) −→ (4, 2) −→ (3, 2) −→ (3, 3) with length 7. However, there might be other permutations, for example, starting at a different interface location, such as (3, 3) −→ (4, 3)−→ (4, 2)−→ (3, 2)−→ (3, 3) consisting of 4 just steps. Each such attempt to construct a permutation σ for a single customer order is called an episode. For each customer order p in historical data P history , we therefore run a specified number of episodes (simulations) that may yield different feasible solutions to the inner problem (5.2). Objective values from solutions across all episodes of each historical order p∈P history are summed and normalized by the number of episodes. In this way, the permutations σ with lower cost have a higher probability of being involved in the calculation of placement quality, while the contribution of less probable but more expensive ones serves as a regularization to derive more robust placement, reflecting the occasional need for "suboptimal" routes at the operational level. 13 The pseudocode for the evaluation of a single placement solution φ is given in Algorithm 1. In compu- tations where the following location to move to (line 9) is sampled, if δ j = 0, then we use the convention that 1/δ j = 1. This essentially states that if on some location we can dispense more than one drug for the considered order p, then all of these drugs are dispensed there (in an arbitrary order), taking advantage of the packing of dispensers that reflects a positive correlation between drug occurrences (see Appendix B). Algorithm 1 Approximate solution to inner SHPP problem: fitness computation. Require: layout H, placement φ :D∪I →H, number of #episodes. Ensure: placement(φ). 1: total← 0 2: for all p∈P history do▷ For each order in historical data. 3: for all e∈1,..., #episodes do▷ Simulate episodes. 4:sample i ⋆ ∼I uniformly▷ Sample uniformly an initial interface tile. 5:R←g|∀(g, ∆ g p )∈ p▷ Drugs requested for order p. 6:repeat 7:L← S g∈R D −1 (g)▷ Set of all dispensers supplying remaining g ∈ R drugs. 8:compute distance δ j =∥φ(i ⋆ )− φ(j)∥ 1 , ∀j ∈ L 9:sample j ⋆ ∼ L proportionally to 1/δ j 10:total← total + δ j ⋆ ▷ Mover travels from location φ(i ⋆ ) to φ(j ⋆ ). 11:R← R\φ −1 (j ⋆ )▷ Remove drug that was dispensed at tile j ⋆ . 12:i ⋆ ← j ⋆ 13:until R =∅▷ Until all drugs have been dispensed. 14:compute distance δ ′ j =∥φ(i ⋆ )− φ(j)∥ 1 , ∀j ∈I ▷ Distance from φ(i ⋆ ) to all interface locations. 15:sample j ⋆ ∼I proportionally to 1/δ ′ j 16:total← total + δ ′ j ⋆ ▷ Go to a finishing interface tile. 17: end for 18: end for 19: return total/ (|P history |· #episodes) ▷ Return the expected number of mover steps in placement φ. This algorithm is used for the computation of the fitness function that steers the metaheuristic algorithm towards placement φ that minimizes (5.1). For this, we use a standard genetic algorithm (GA) that represents placements φ : D∪I → H as permutations that map (indexed) packed tiles D to specific positions in the layoutH. The algorithm is implemented in the Pymoo library (Blank & Deb, 2020) and uses order crossover with inversion mutation operators. Furthermore, we would like to highlight the fact that the Algorithm 1 relaxes on the dispensing time; thus, it acts as if it spent zero time at each dispensing site. Hence, congestion at the dispensers is not taken into account by the placement—all orders are essentially being completed simultaneously. However, this is handled indirectly by sampling next dispensers in the sequence σ to visit proportionally to the inverse distance from the current location. This aims to reflect that the nearest compatible dispenser may not always be available. Finally, we also tested a different placement algorithm, known as the Analytical placer, which solved the inner optimization problem (5.2) optimally by enumerating all permutations σ. This is computationally tractable if the number of alternative dispensers z g for all drugs g ∈G is not too high, as the typical order consists of approximately 4 drugs. However, with more alternatives in dispensers (larger neighborhoods), this quickly leads to prohibitive computational time. What is perhaps more important, however, is the fact that the placement objective (5.1) is sensitive to the best-case traversal through the given layout due to 14 the inner problem (5.2). This means that, e.g., the second-shortest routes are not optimized, although they might be frequently used in scheduling at the operational level, especially when more than one mover is used. 6. Scheduling problem The Scheduling problem addresses the operational-level decisions of executing a real, given set of ordersP on the given layoutH with a placement φ of dispensersD and interfacesI. The goal is to construct a schedule s of dispensing operations for the given set of orders P. The schedules assigns all dispensing operations of orders to movers while adhering to constraints: the capacity of the movers must not be exceeded, and each dispenser and tile can be occupied by only one mover during dispensing operation. The mover realizes one order p∈P at a time, starting and finishing at some interface locationI. All operations must be realized at the physical locations of the corresponding dispensers and interfaces as defined by the placement φ. Travel time between two operations performed at different tiles is equal to their ℓ 1 distance in the layout H. In scheduling terms, each order p ∈ P represents a job. Every order p consists of operations T p . Each operation is either dispensing of drug g (with duration ∆ g p ) or a cartridge exchange (loading/collection - both with constant duration η I ). Interfaces and drug dispensers are represented as resources; every dispenser for drug g can realize any g-dispensing operation, and every interface can realize any cartridge exchange operation. Two more resources are consumed with each operation - the tile below the assigned dispenser and a mover m. Mover m is a take-give resource, meaning all operations from an order p must be executed successively on the same mover m before another order p ′ can be realized by m. Travel times are incorporated through a transition distance matrix (later defined as L tiles in Section 6.2.2) that specifies the minimum time interval required between consecutive operations based on their relative distance. There are no precedence constraints between the operations within an order, except that cartridge loading must be the first operation and cartridge collection the last. The output of the scheduling algorithm comprises the assignment of each order to a specific mover, the start times of all operations, the allocation of each drug-dispensing operation to a particular dispenser (in cases where multiple dispensers are available), and the assignment of each start and finish (cartridge loading and collecting, respectively) operation to specific interfaces. Formally, we define a set of all individual dispensing operations and cartridge exchange operations T = S p∈P T p , where for each order p ∈ P we define a set of 4-tuples, one for each drug g ∈ S p in order p, one for a start operation and one for a finish operation: T p =(p,g, ∆ g p , dispensing)|∀g ∈S p ∪(p, interface,η I , start), (p, interface,η I , finish). Thus, each operation is defined by a 4-tuple (p,g, ∆ g p ,τ ), where p∈P denotes the order, g ∈G∪interface is the drug or interface, and ∆ g p ∈N is the duration (interface operations always require η I time units). The type τ distinguishes between dispensing drug g ∈G and the loading (start) or collecting (finish) of cartridges. 6.1. Constraint Programming Model Variables Generally, scheduling problems formulated as Constraint Programming models use the concept of interval variables (Heinz et al., 2022), which serve as the primary building blocks of the model. For each 4-tuple t∈T , 15 we create an interval variable task_var t with a duration ∆ g p or η I representing its execution, and whose start time is determined by the scheduler. To allow the solver to decide which dispensing operation will be performed on which mover and at which dispenser, we define the set of alternatives as A = S t∈T A t , where for each operation t = (p,g, ∆,τ ) ∈ T , we define the set A t = (m,t,h) | ∀m ∈ M, ∀h ∈ H : g ∈ φ −1 (h), where φ −1 (h) is the set of drugs dispensed at tile h∈H. Each alternative (m,t,h) corresponds to the execution of operation t by mover m at tile h, while pre- serving the duration ∆ g p and type τ of the original operation. We then create interval variables for all alternative tuples in the same manner as with the original tuples t∈T . We denote these interval variables alt_task m,t,h . To ensure that two interval variables do not overlap, they must be included within the same wrapper (data structure known as a sequence variable in IBM CP Optimizer). Consequently, multiple wrappers are used to group interval variables based on shared properties. An illustration of these wrappers is provided in Figure 9. It shows how interval variables for a specific order on a specific mover are grouped (e.g., order_wrap m 1 ,p 6 ), and how these order-specific groups are then collected under a main wrapper for each mover (e.g., mover_wrap m 0 ). VariableDescription task_var t Interval variable, ∀t ∈ T . Represents filling of a specific drug for an order or loading/unloading its cartridge. alt_task m,t,h Optional interval variable, ∀(m, t, h) ∈ A. Represents the possibility of mover m executing operation t on tile h. order_wrap m,p Optional wrapper, ∀m ∈ M, ∀p ∈ P. Represents order p being served by mover m. Spans across all interval variables for order p. wrapper_mover_tasks m Wrapper for all operations assignable to mover m. For example, for m ′ ∈ M, wrapper_mover_tasks m ′ = alt_task m,t,h | ∀(m, t, h) ∈ A, m = m ′ . wrapper_identical_tiles h Wrapper for all operations executable to tile h. For example, for h ′ ∈ H, wrapper_identical_tiles h ′ = alt_task m,t,h | ∀(m, t, h) ∈ A, h = h ′ . mover_wrap m Wrapper: mover_wrap m,p | ∀p ∈ P, ∀m ∈ M. Wrapper for all orders served by mover m. Table 1: Scheduling model variables. 6.2. Constraint Programming Model Constraints The optimal allocation of dispensing operations t ∈ T to movers lies at the core of this scheduling problem. The objective is to assign operations to movers so that the total execution time is minimal while adhering to several vital constraints. 6.2.1. Alternative Movers and Dispensers For every operation t∈T , we define task_var t as a master interval variable. For each feasible (mover, tile) execution option, we have defined optional alternatives alt_task m,t,h . For every t ′ ∈ T , we link the master variable task_var t ′ to its alternatives using alternative (task_var t ′ ,alt_task m,t,h |∀(m,t,h)∈A,t = t ′ ).(6.1) This guarantees that exactly one alternative is chosen for the given task_var t ′ . 16 mover_wrap m 0 mover_wrap m 1 order_wrap m 0 ,p 0 order_wrap m 0 ,p 5 order_wrap m 1 ,p 6 Input: Pool of Alternatives (alt_task) ∈ A Scheduler S hctz simva F S lisino F ... S simva lisino F ... wrapper_identical_tiles 1×1 Output: Optimized Schedule Figure 9: The scheduler maps input operations (left) to an optimized assigned schedule (right). Depicted intervals represent chosen operation alternatives for movers m 0 and m 1 . This schedule fulfills three orders via order_wrap m 0 ,p 0 , order_wrap m 0 ,p 5 , and order_wrap m 1 ,p 6 . Both lisinopril operations are grouped under wrapper_identical_tiles 1×1 , indicating both were dispensed from the same dispenser at tile (1,1). The two simvastatin operations are not grouped in such a way, indicating each was filled at a different dispenser. 6.2.2. No Overlap of Mover Activities Our initial step in model construction is to prohibit movers from using multiple tiles at the same time, i.e., a single mover cannot be in multiple physical locations at once. We enforce no overlap among all interval variables targeting mover m ∈ M. Furthermore, we enforce travel times (i.e., gaps in the schedule) between any two interval variables alt_task m,t,h and alt_task m ′ ,t ′ ,h ′ , ∀(m,t,h), (m ′ ,t ′ ,h ′ ) ∈ A, equal to the physical distance between the corresponding dispensers’ positions. The distance matrix L tiles captures the travel time between tiles h,h ′ ∈H, defined as: L tiles [(m,t,h), (m ′ ,t ′ ,h ′ )] =∥h− h ′ ∥ 1 .(6.2) We pass L tiles to the model as part of no_overlap constraint to enforce the travel times: no_overlap (alt_task m,t,h |∀(m,t,h)∈A,m = m ′ ,L tiles ) ∀m ′ ∈M.(6.3) 6.2.3. Order Completion Ensuring non-overlapping execution of dispensing operations is a necessary condition but not sufficient. We also have to enforce that no two dispensing operations of different orders interleave; an order must be completed before starting a new one. For each order p ′ ∈P and for each mover m ′ ∈M, we group the corresponding optional interval variables into a wrapper order_wrap p ′ ,m ′ , defined as: order_wrap p ′ ,m ′ =alt_task m,t,h | (m,t,h)∈A,m = m ′ ,p ′ ∈ t. Then, for each mover m ′ ∈ M, we group all order wrappers order_wrap p,m ′ for all p ∈ P into a mover_wrap m ′ : mover_wrap m =order_wrap p,m |∀p∈P ∀m∈M. 17 This hierarchical structure is visualized in Figure 9. For instance, mover m 0 is assigned orders p 0 and p 5 , thus the top-level wrapper mover_wrap m 0 contains the individual order wrappers order_wrap m 0 ,p 0 and order_wrap m 0 ,p 5 . As seen in Figure 9, each order_wrap wrapper begins with a start operation and concludes with a finish operation. To ensure a feasible ordering, we must guarantee that all dispensing interval variables occur strictly between the start and finish intervals. To formalize this, we first define three subsets of alternative tasks for each order p∈P and mover m∈M. Let Alt start p,m be the set of operations representing loading and empty cartridge,Alt finish p,m be the set of operations representing unloading a completed cartridge, and let Alt dispensing p,m be the set of operations representing drug dispensing: Alt τ p,m =alt_task m,t,h | (m,t,h)∈A, p∈ t∧ τ ∈ t ∀τ ∈start, dispensing, finish. We then apply a constraint that ensures every possible start operation begins before any dispensing oper- ation. This guarantees that regardless of which interface is selected, the initialization always precedes any dispensing: start_before_start (a start ,a disp ), ∀p∈P,∀m∈M,∀a start ∈Alt start p,m ,∀a disp ∈Alt dispensing p,m (6.4) Similarly, we must enforce that every dispensing operation ends before the finish operation ends. This prevents any drug filling task from occurring after the order wrapper has supposedly finished: end_before_end (a disp ,a finish ), ∀p∈P,∀m∈M,∀a disp ∈Alt dispensing p,m ,∀a finish ∈Alt finish p,m (6.5) This yields a result comparable to what we would obtain using a span constraint. Finally, we make sure that dispensing operations for different orders assigned to the same mover do not overlap: no_overlap (mover_wrap m ) ∀m∈M.(6.6) This not only enforces non-overlapping execution, but also guarantees full separation between orders: a mover must complete all operations for an order p before beginning any loading/dispensing operations for another order p ′ . Additional mechanisms are required to ensure that the previous constraint operates as intended. In its current formulation, each mover wrapper contains order wrappers corresponding to all orders. The following technique ensures that an order wrapper is instantiated only when the corresponding mover–order assignment alternative is selected (i.e., if order 3 is assigned to mover 2, then only the wrapper order_wrap 3,2 is instantiated and placed within mover_wrap 2 ). For each order p∈P, we define u t ′ ,m ′ = X a∈B presence(a) ∀t ′ ∈T p ,∀m ′ ∈M,B =alt_task t,m,h | (t,m,h)∈A,t = t ′ ,m = m ′ . For each operation t ∈ T p and mover m ∈ M, the variable u t,m reflects whether mover m executes operation t: it equals 1 only for the mover assigned to that operation and 0 otherwise. Let t 1 be the first operation in T p . By enforcing u p,m = u t 1 ,m we propagate the existing mover assignment 18 information from the operation-level variables u t,m to a single representative variable u p,m for the entire order p — u p,m will be equal to 1 only for the mover executing the order and 0 for all others. We then make sure order_wrap p,m is present only when u p,m is 1: presence(order_wrap p,m ) = u p,m . 6.2.4. No Overlap over Dispensing Activities Every dispenser must also allow at most one operation at any given moment. To achieve that, for any given tile h ′ ∈ H, we create a wrapper sequence_identical_tiles h ′ of alt_tasks m,t,h ,∀m ∈ M,∀t ∈ T ,h = h ′ . Figure 9 shows an example of this with sequence_identical_tiles 1×1 . This wrapper groups interval variables across different orders (p 5 ,p 6 ) and movers (m 0 ,m 1 ) that all require the same dispenser tile. We ensure the tile is not used concurrently by multiple movers with another no_overlap constraint: no_overlap (wrapper_identical_tiles h ) ∀h∈H.(6.7) Finally, the model’s objective is to minimize the latest finish time across all interval variables: min maxend_of(task_var t )|∀t∈T.(6.8) This completes the Constraint Programming formulation for the core scheduling problem, which assigns dispensing operations to movers and determines their order. The model captures all hard constraints— mover, tile, and dispenser exclusivity, travel times between dispensers, while intentionally relaxing conflict resolution among movers sharing pathways. In the following subsection, we present a calculation of a lower bound on Scheduling. 6.3. Lower Bound: Parallel Movers Using SHPPN We propose a lower bound on the optimal makespan with the two goals: (a) to estimate the quality of the solution of Scheduling, and (b) to provide a warm start solution for the solver. The main idea of the lower bound consists of relaxing the overlaps on dispensers and the travel times between interfaces. On the other hand, we include travel times between the dispensers themselves. This is achieved by constructing an instance of a TSP with neighborhoods (TSPN) for each order p∈P. Subsequently, the results from TSPN would be used as constants in the parallel machine scheduling problem P||C max , with jobs having processing times set to the lengths of the optimal tours found by TSPN and a set of resources corresponding to the available movers. Let κ p =min σ∈I×Π(p)×I |σ|−1 X j=1 ∥φ(σ j )− φ(σ j+1 )∥ 1 (6.9) be the length of the shortest path between some two interface locations that visits all dispensersS p required for order p. Then we define η ′ p = 2η I + κ p + X (g,∆ p g )∈p ∆ p g ,(6.10) denoting the lower bound on the total processing time for order p, where η I is the time of releasing/col- lecting an empty cartridge, κ p is the shortest possible traveling time of the mover performing order p, and ∆ p g is the time to dispense drug g in customer order p. 19 The values η ′ p can be computed by solving one instance of the SHPPN for each customer order p sepa- rately. However, SHPPN can be reduced to the TSPN problem and subsequently reduced to standard TSP using the reduction of Noon & Bean (1993) with the introduction of only a single additional vertex. Since the instance of the SHPPN typically involves only small graphs (e.g., ≈ 30 vertices since a typical order might consist of 4 drugs plus additional start and finish interface locations, each with approximately 5 alter- natives), the resulting instances of standard TSP can be solved efficiently by, e.g., lazy subtour elimination in a fraction of a second. With the values of η ′ p , the lower bound is computed as follows. Let us define an instance of the P||C max problem (i.e., scheduling on parallel identical machines with makespan minimization) with n movers resources and a set of n jobs with processing times η ′ 1 ,...,η ′ |P| . We claim that P,n movers |η ′ p n i=1 |C max acts as a lower bound on schedule makespan. Proposition 2 (SHPPN lower bound). Optimal solution of P,n movers |η ′ p n i=1 |C max is a lower bound on Scheduling. Without detailing the proof, the main idea of SHPPN lower bound is that the optimal solution of P,n movers |η ′ p n i=1 |C max satisfies all constraints on Scheduling, detailed in Section 6.2, except it relaxes on: (i) the overlap at dispensers, tiles and interfaces locations, and (i) the travel time between two interface locations where a mover finishes with a completed order and the following interface where it starts with the next order (i.e., instantaneous teleportation between interface locations). Regarding the complexity of computing the SHPPN lower bound, it is known that P||C max is (weakly) NP-hard already for two machines (i.e., movers) and strongly NP-hard for an unbounded number of machines. Computing values of η ′ p for the realistic scale of the problem can be done efficiently in practice with an existing TSP solver. Therefore, we deem that the lower bound is computationally practical. Additionally, we also use the solution computed by the lower bound procedure to warm-start the Constraint Programming model for the Scheduling problem by supplying mover-order assignments, which accelerate its solution. 7. Routing Problem The objective of the Routing problem is to construct feasible ℓ 1 paths for each mover m∈M based on the pre-determined schedule of dispensing operations. Our routing focuses on resolving conflicts, defined here as simultaneous tile occupancy. Unlike classical MAPF, our routing has to precisely follow the schedule. Furthermore, spatial overlap of two movers on one tile is permissible provided no dispensing occurs. If during dispensing on one mover another has to travel through the same tile, the dispensing is paused. Since standard MAPF solvers cannot accommodate these context-dependent rules, we employ a tailored iterative heuristic. The algorithm begins by generating fixed resting sites — locations where movers can wait without obstructing access to dispensers (Section 7.1). In each iteration, an ILP assigns idle movers to these sites (Section 7.2), and a polynomial-time algorithm constructs least-congested shortest paths between operations and resting spots. These paths are used to detect conflicts. Although rare, if any arise, the schedule is revised and the process repeats if needed. 20 7.1. Resting Sites Generation Movers may experience idle times between operation executions due to dispenser unavailability. To manage these times efficiently and maintain system throughput, we assign resting sites — specific locations where movers can remain stationary without interfering with dispensing operations on adjacent tiles. We consider a fixed dispenser nozzle at the center of each tile, which enables cooperative operation by placing resting sites at tile intersections. This layout allows a resting mover (represented in gray in Figure 10) to remain idle while another mover utilizes an adjacent tile for transit or dispensing. To maintain clearance for these resting movers at the intersections, the strategy leverages the mover’s rotational capabilities: after filling the capsules on one side of the cartridge (purple side, Figure 10a), the mover shifts slightly to rotate 180 ◦ (depicted in Figure 10b) and fill the other side (green side, Figure 10c), which together comprise the full cartridge. By utilizing this rotation instead of traversing the entire tile to fill every capsule, the mover minimizes its footprint (the patterned area in Figure 10), thereby preserving the necessary clearance for half a resting site at the boundary. Resting MoverLoading Mover (a) Purple Area Can Be Filled Resting Mover Loading Mover (b) 180° Rotation Resting Mover Loading Mover (c) Green Area Can Be Filled Figure 10: Visualization of dispensing maneuver. Geometrically, a mover with dimensions a×a rotates within a disk of radius R = a/ √ 2. On a 2a× 2a tile, this leaves free space of 2a− 2R = a(2− √ 2) > 0.5a along the edge, which is sufficient to accommodate half of a resting mover. However, the rotation depicted in Figure 10 imposes a constraint: each dispensing tile supports only one resting site to preserve rotation space for the dispensing mover. Interface tiles are an exception, accommo- dating two resting sites since movers remain centered during cartridge swaps. We formulate the problem of finding viable resting sites layout as a Maximum Independent Set problem over a graph (see Figure 11b) constructed from layoutH. In this graph, edges (depicted as red line segments) represent mutual exclusion constraints — that is, two nodes (potential resting sites, depicted as blue dots) are connected if resting at both simultaneously would cause dispensing inoperability on their common tile. The resulting solution, shown in Figure 11c, highlights the selected resting sites R using yellow rectangles. (a) Layout H with resting sites candi- dates in blue. (b) Constructed Maximum Independent Set graph. (c) Solution for Maximum Independent Set on a 4× 4∼ 2 square layout. Figure 11: The yellow areas represent selected resting sites R for the movers out of the candidate positions depicted by the blue dots. Interface tiles (gray, hatched tiles) may contain up to 2 resting sites. 21 AMLODIPINEATENOLOLATORVASTATINFUROSEMIDEHYDROCHLOROTHIAZIDE LEVOTHYROXINELISINOPRILLOSARTANMETFORMINOMEPRAZOLE SIMVASTATINVALSARTANINTERFACE Mover 0 Mover 1 Mover 2 Mover 3 0100200300400500600700800 (a) Optimized schedule colored by drug type. Mover 0 Mover 1 Mover 2 Mover 3 0100200300400500600700800 w 1 w 2 w 3 w 4 (b) Transits with idle periods - input for the routing phase. Mover 0 Mover 1 Mover 2 Mover 3 0100200300400500600700800 r 1 r 2 r 6 r 1 (c) Transits with idle periods with assigned resting sites. Figure 12: Optimized schedules, identification of its idle periods and their assignment 7.2. Resting Sites Assignment The problem of assigning each mover’s idle period to a resting site in the schedules is formulated as a Mixed-Integer Programming model. The model takes as input the set of precomputed resting sites R = r 1 ,r 2 ,...,r k , the set of moversM =m 1 ,m 2 ,...,m n movers , and the set of transits W =w 1 ,w 2 ,...,w l derived from the schedules. A transit represents the time gap between two consecutive operations assigned to the same mover in the schedule. Among these, transits with idle periods are those whose duration exceeds the travel time required to move between the corresponding operation locations (see Figure 12b). We populate set W only with transits with idle periods, as these represent the intervals during which a mover needs a resting site. The objective of the model is to assign each such transit to a resting site in a way that minimizes the additional travel cost for movers to reach and leave the resting sites. The problem is modeled using binary decision variables x ij ∈ 0, 1, where x ij = 1 if transit w i ∈ W is assigned to resting site r j ∈ R, and c ij denotes the round-trip travel time. The objective is to minimize P i,j c ij x ij while assigning each transit to a resting site: P j∈R x ij = 1, ∀i∈W and ensuring at most one mover uses a resting site at any given time: x ij + x i ′ j ≤ 1, ∀(a ij ,a i ′ j ) ∈ F. Set F consists of pairs of operations with overlapping time windows. 7.3. Conflict Resolution Unlike standard MAPF collisions, our system allows movers to coexist on a tile. A conflict occurs specifically when a mover transits through a tile where another is actively dispensing, forcing the operation 22 IFHCTZSIMVLISINIFIFLISINSIMVATENFUROIF k 0 IFSIMVLISINMETFIFIFMETFOMEPLOSARIFIFLISINMETFOMEPIF IFLISINAMLODFUROLEVOMETFSIMVIF IFMETFATORVHCTZIFIFSIMVLISINVALSAHCTZIF Figure 13: DAG generated from the schedule in Figure 12a. Black edges represent same-mover precedence constraints, blue edges represent same-dispenser precedence constraints. Task colors indicate dispenser assignment. to pause (see Figure 10, loading mover has to yield so there is space for transit). Our algorithm detects these interruptions and extends the schedule to ensure the required dispensing duration is met. We model the schedule as a directed acyclic graph (DAG) Ω = (K∪k 0 ,E), where vertices K represent operations and k 0 is a virtual source node. Edges E represent start-to-start precedence constraints: 1. Same-mover constraint: Connects sequential operations of a single mover. We include edges (k 0 ,k) for the first operation k of each mover with a weight equal to the start time of k. For all other mover edges, the weight w(k,k ′ ) = l i + ∆ g ′ p ′ +∥h ′ − h ′ ∥ 1 accounts for the interruption time l i , dispensing duration ∆ g ′ p ′ , and travel time between locations h ′ and h ′ . 2. Same-dispenser constraint: Connects operations sharing a dispenser, ordered by their original start times, with a weight of 1. To adjust the schedule, we calculate the longest path from the source vertex to every vertex. Setting the start time of the source node S(k 0 ) = 0, the start time S(k) for all other operations is determined recursively: S(k) = max (k ′ ,k)∈E S(k ′ ) + w(k ′ ,k). The global makespan is then determined by max k∈K S(k) + ∆ g ′ p ′ . This process is iterative. We recalculate the interruptions and update the new interruption counts using l i = maxl i , current interruptions and re-optimize the start times if any values have changed. If not, we have reached a conflict-free solution. 7.4. DAG-based Schedule Merging For larger order sets, it might be useful to split the set into smaller batches randomly and independently generate sub-schedules for each batch. These sub-schedules are then transformed into DAGs using the same principle and connected together. Specifically, edges are added from the terminal vertices of the preceding DAG to the initial vertices of the new sub-schedule. Edge weights correspond to the minimum required transit times between the last and first tasks of adjacent batches. Additionally, edges are added between the last tasks in the earlier schedule and the first tasks in the new schedule that use the same tile, ensuring potential dispensing conflicts are accounted for. The start times of the new intervals are then initialized by propagating earliest start times through the extended graph, computed via LP or longest-path relaxation over the merged DAG. 23 8. Experiments The proposed algorithms allow investors to evaluate potential facility configurations — including tile layouts, dispenser placements, and mover counts — and assess whether a given investment yields sufficient production throughput. This section pursues two aims: validating the proposed methodology and demon- strating that the decomposition introduces no significant loss in solution quality relative to the theoretical lower bound (Section 8.6). At the tactical level, we evaluate robustness (Section 8.2) and analyze how tacti- cal decisions propagate to the operational level (Section 8.3). At the operational level, we examine the effect of mover count and time limits on makespan (Sections 8.4 and 8.5), assess the scalability of DAG-based schedule merging and its effect on solution quality for larger problem instances (Section 8.7), and quantify the impact of the routing phase on overall system performance (Section 8.8). As a basis for our study, we applied data from an analysis of real-world co-prescription patterns within the NHANES drug survey dataset (National Center for Health Statistics, 2012). The analysis provided the marginal probabilities and correlation matrix (discussed in Packing, illustrated in Figure B.1 for the top 40 most frequently prescribed drugs. To reflect sufficient system complexity, simulated orders were restricted to contain between 3 and 8 unique medications. To generate realistic, typical order-sets, we implemented a synthetic data generator that produces mul- tivariate binary patient prescriptions. This system is based on the methodology proposed by Leisch et al. (1998), which transforms samples drawn from a multivariate normal distribution into binary outcomes via marginal thresholding. The generated instances can be found on GitHub [link removed while under review]. 8.1. Hardware and Software Environment All computational experiments were conducted on a server equipped with dual Intel Xeon Silver 4110 CPUs @ 2.10GHz, providing a total of 16 physical cores and 32 logical threads. The system is provisioned with 187 GB of RAM. The operating environment is Debian GNU/Linux 13 (trixie). The algorithms were implemented in Python 3. Mathematical optimization models were solved using commercial solvers, with the environment configured to support Gurobi 12.0.0 and CPLEX 22.1.1.0. To ensure reproducibility, we provide a summary of the parameters used for our algorithms in Table 2. The dispensing speed parameter we derived from Sonntag et al. (2023). As future works are primarily aimed at speeding up the dispensing process, we also evaluated the system’s performance with faster dispenser speeds. Tactical Level (Placement & Packing)Operational Level (Scheduling, Routing) GA Population Size 150Time Limit (Standard)600s GA Max Evaluations 50,000Time Limit (Scalability)300s – 2,000s GA Episodes20Standard Dispensing Speed 100 ticks / cartridge Dispenser Count82 (Fig. 14), 62 (Fig. 15) Fast Dispensing Speed10 ticks / cartridge Table 2: Summary of experimental parameters. 24 8.2. Robustness of Tactical Level To ensure that the placement and packing phases do not suffer from overfitting, where the resulting layout is optimized for a specific sequence of orders and, therefore, performs poorly on others, we evaluated the robustness of the solution. The packing and placement phases were run using a training set of 100 generated orders to create a fixed 8× 8 layout. The packing configuration consisted of 82 dispensers, optimized by a genetic algorithm with a population size of 150 and a maximum of 50,000 evaluations. Subsequently, to test the robustness of the resulting layout, we executed the scheduling and routing phases using 100 distinct datasets: the original training set and 99 new, unseen validation sets, each consisting of 100 orders. For these validation runs, the scheduler used a warm-start initialization and was constrained to a 600-second time limit. The results in Figure 14 demonstrate a high degree of stability across the different datasets. The −10−50510 0 5 10 15 20 relative difference [%] count [-] (a) Relative makespan deviation from the mean across 100 runs, for a layout with 8 movers. −10−50510 0 5 10 15 20 relative difference [%] (b) Relative makespan deviation from the mean across 100 runs, for a layout with 12 movers. 0 10 0 10 y x (c) Interface occurrence heatmap with 2 interfaces in each layout type. From top to bottom, the lay- outs doubleline, line, ring, and square are displayed. Figure 14: Robustness and stability analysis of the tactical level across varying mover configurations and layout topologies. makespan (C max ) observed for the unseen validation sets deviates by no more than ±10% from the result obtained with the training set. This narrow margin indicates that the layout generated by the placement and packing phases is robust and not over-fitted to the particular historic set P history . As a further experiment to assess the robustness of the tactical level, the stability of the interface placement is assessed. The conducted experiment consisted of 100 independent placement optimizations with generated sets mirroring the medicine distribution of historic order P history on line, ring, doubleline, and square layouts to identify the most occurring interface locations when 2 interfaces are used. A stable placement should put interface tiles into similar positions when repeated with similar inputs. The results in Figure 14c show the interface heatmap for all 4 considered layouts, with color indicating the frequency of interface placement. Due to the rotation symmetry of the ring topology, the resulting layouts were canonized so that their first interface is placed at the same coordinates. The results suggest a strong preference for clustering interfaces at the geometric center of the grid, while the peripheral tiles were used almost exclusively for dispensers. While this configuration minimizes theoretical mover travel times, it may pose practical engineering challenges regarding the physical accessibility of the interface for cartridge loading and unloading. 25 8.3. Correlation of Placement and Scheduling Objectives To validate the link between the proposed objective function for Placement and Scheduling, we analyzed the correlation between the objective values (5.1) obtained during placement optimization via GA and the resulting final scheduling makespan (C max ). The experimental setup used layouts with 48 tiles (47 for the square topology) and 2 interfaces, with 5 movers fulfilling 50 orders. The scheduling phase was conducted without warm-start initialization, constrained by a 600-second time limit across 5 CPU threads, with a packing configuration of 62 dispensers on the fast dispensing setting (see Table 2). 101520253035 400 450 500 550 600 placement objective [-] C max [-] Line Doubleline Ring Square Figure 15: Correlation of placement and scheduling objective. The results, illustrated in Figure 15, demonstrate the progression of the GA across different topologies. Each data point represents a placement solution evaluated during the Placement optimization. The scatter plot reveals a linear relationship: as the GA minimizes the placement objective (5.1) through successive generations, the corresponding C max decreases proportionally. This consistent correlation confirms that the placement objective (5.1) serves as a reasonable proxy for the actual scheduling. By optimizing for this proxy, Placement effectively drives the search toward configurations that minimize the final execution time. Furthermore, the model correctly identifies the performance potential of different structural layouts, accurately placing the square configuration in the highest-efficiency region (lowest C max ) and simultaneously highlighting the inherent physical constraints and higher costs associated with the ring and line topologies. 8.4. Effect of Mover Count on C max Figure 16 illustrates the makespan trend as the number of movers increases across the four considered layouts. Initially, adding more movers significantly reduces the makespan, since operations can be distributed more evenly across the system. However, the improvement rate diminishes as the number of movers grows. This plateau effect is explained by the limited number of available dispensers per drug type; adding additional movers beyond a certain point does not yield further speed-up, as the dispensers themselves become the throughput bottleneck. The exact makespan values for each configuration are reported in Tables 4, 5. 8.5. Effect of Time Limit on C max The scheduling problem grows exponentially in complexity with the number of orders (if the batching mechanism detailed in Subsection 7.4 is not used), making exact optimization infeasible for even moderately sized instances. The aim is therefore to obtain high-quality solutions within practical time limits, raising the question of when to terminate the search. Figure 17a evaluates this trade-off by running the scheduler on order sets of sizes 10, 25, 100, and 250 under varying time limits, while Figure 17b extends the analysis to sets of 200, 300, and 500 orders. 26 123456789101112 0 10,000 20,000 30,000 number of movers [-] C max [-] Line Ring Doubleline Square (a) Standard dispensing speed. 123456789101112 0 2,000 4,000 6,000 number of movers [-] Line Ring Doubleline Square (b) Fast dispensing speed. Figure 16: Effect of mover count on C max for different topologies (64 tiles each) with varying dispensing speeds. 01,0002,000 0 2,000 4,000 6,000 time [s] C max [-] 01,0002,000 0 10,000 20,000 30,000 time [s] |P| :10 2050 100200 300500 Figure 17: The effect of the time limit on C max for different sizes of P on a square 8× 8∼ 2 topology. The results for order sets with 100 or fewer orders show diminishing improvements in makespan with longer run times. Time limits of around 300 seconds appear to be sufficient for practical reasons across all order set sizes. Figure 17b illustrates the impact of instance size on problem complexity. For larger order sets, even identifying a feasible solution requires significantly more time, and subsequent improvements become neg- ligible within any reasonable time limit. This behavior stems from the exponential growth of the solution space that the solver must explore. 8.6. Comparison Between C max and Lower Bound Because the previous analysis focused on identifying practical time limits by comparing makespan across different solver runs, it only provided a relative notion of solution quality. To put these results into a more absolute context, we would ideally compare the obtained schedules against the theoretical optimal makespan. Since computing the exact optimum is computationally intractable, we instead use the lower bound introduced in Section 6.3. As this bound is derived by relaxing certain problem constraints, the true optimum must lie strictly above it. Consequently, the reported gaps in Table 3 are conservative, and the actual distance of our solutions from optimality is even smaller than the numbers suggest. The results show that the C max of our solutions is consistently close to the lower bound, with gaps generally remaining in a moderate range. For larger order sets combined with many movers, the differences increase more noticeably, reinforcing that problem complexity at this scale begins to limit how close the solver can remain to the optimum. 27 Mover Count Order Set Size: 25Order Set Size: 50Order Set Size: 100 LB C max Gap Gap [%]LB C max Gap Gap [%]LB C max Gap Gap [%] 23804 4173 369 8.84%7541 8518 977 11.47%15304 16727 1423 8.51% 6 1270 1577 307 19.47%2514 3008 494 16.42%5101 5879 778 13.23% 10767 1010 243 24.06%1509 1896 387 20.41%3063 3692 629 17.04% 12–2555 3416 861 25.20% Table 3: Comparison of lower bound vs. actual makespan across different order set sizes (8× 8∼ 2 layout). 8.7. Achieved Scalability: Effect of DAG-based Schedule Merging Subsection 7.4 discusses batching larger order sets and then merging them together. To measure the impact this technique has on the overall C max , we test different batch sizes for varying numbers of orders while keeping the number of movers fixed at eight. Each curve shows the evolution of the makespan (C max ) over time as the solver continues to improve the solution. 01,0002,000 4,200 4,400 4,600 time [s] C max (a) |P| = 100. 01,0002,000 8,500 9,000 9,500 time [s] (b) |P| = 200. 01,0002,000 13,000 14,000 15,000 time [s] C max (c) |P| = 300. 01,0002,000 22,000 24,000 26,000 time [s] (d) |P| = 500. Batch size: 50 100 200 300 500 Figure 18: The effect of the time limit on C max for different batch sizes. Square 8× 8∼ 2 layout, 8 movers. The results demonstrate that partitioning the problem into smaller batches is an effective strategy for scaling to larger order sets. Although this decomposition introduces minor suboptimalities for smaller instances, it substantially improves computational efficiency. For larger instances, the method attains so- lutions of comparable or superior quality (Figure 18d, stitching batches of size 200 (red line) has better performance than no batching (purple line)) within a significantly shorter time horizon, indicating that the trade-off between optimality and tractability is highly favorable under realistic time constraints. 8.8. Effect of the Routing Phase on Solution Quality To assess the effect of the routing phase on system performance, we report the makespan values obtained from the scheduling model and after routing with conflict resolution. Tables 4, 5 present the results for four layout types (ring, square, doubleline and line) across different numbers of movers and two different drug 28 dispensing speeds. The values illustrate the relative changes in makespan between the scheduled and routed solutions. The results indicate that the routing overhead is minimal across all layout types and numbers Mover count Square C max Doubleline C max Ring C max Line C max SchedulingRoutingSchedulingRoutingSchedulingRoutingSchedulingRouting 133427+0 (0.0%)34506+0 (0.0%)35056+0 (0.0%)35233+0 (0.0%) 48373+10 (0.1%)8638+21 (0.2%)8806+29 (0.3%)8875+38 (0.4%) 84224+24 (0.6%)4399+37 (0.8%)4512+46 (1.0%)4560+53 (1.2%) 12 2942+23 (0.8%)3042+42 (1.4%)3108+56 (1.8%)3286+100 (3.0%) Table 4: Makespan values for different movers, before and after routing; standard speed. Mover count Square C max Doubleline C max Ring C max Line C max SchedulingRoutingSchedulingRoutingSchedulingRoutingSchedulingRouting 14731+0 (0.0%)5731+0 (0.0%)6202+0 (0.0%)6201+0 (0.0%) 4 1199+10 (0.8%)1492+23 (1.5%)1714+42 (2.5%)1638+50 (3.1%) 8 619+23 (3.7%)858+39 (4.5%)914+52 (5.7%)980+53 (5.4%) 12448+37 (8.3%)605+56 (9.3%)665+60 (9.0%)689+79 (11.5%) Table 5: Makespan values for different movers, before and after routing; fast dispensing speed. of movers. This demonstrates that the majority of the routing problem is already effectively handled by the scheduling algorithm, confirming that our decomposition approach is well-structured. Table 5 demonstrates that increasing the dispensing speed leads to a larger routing overhead; as the movers spend comparatively more time in transit than dispensing, the likelihood of spatial conflicts increases, requiring more resolution time. 9. Conclusion This work addresses the shift toward mass personalization in pharmaceutical manufacturing by framing the production process as a complex logistics challenge. Inspired by planar transport technologies like the Beckhoff XPlanar, we designed a framework that splits decision-making into a tactical layer for long-term dispenser layout and an operational layer for daily execution. At the tactical level, we optimized the physical configuration of the production line by utilizing a Mixed- Integer Quadratic Programming model to solve the Packing problem, which allocates dispensers based on drug co-occurrence patterns in historical patient data. Subsequently, a Genetic Algorithm addressed the Placement problem by arranging these dispensers to minimize expected mover travel distances. Our experiments demonstrated that the objective function used in this phase strongly correlates with the final scheduling makespan, while robustness analysis confirmed that these layouts perform consistently well on unseen validation sets. At the operational level, we focused on daily order fulfillment by formulating the Scheduling problem through Constraint Programming, modeling movers as take-give reservoir resources. To ensure scalability, we handle larger order sets via batch decomposition. The Scheduling phase is complemented by a Routing phase that uses iterative conflict-resolution and DAG-based logic to generate collision-free paths. Results indicate that routing overhead remains minimal, confirming that the schedule accurately approximates physical 29 movements. Practically, this framework serves as a simulation tool for investors and system designers to estimate throughput and identify diminishing returns on hardware resources before physical implementation. Declaration of generative AI and AI-assisted technologies in the manuscript preparation pro- cess. During the preparation of this work, the authors used ChatGPT and Gemini to assist with refining linguistic clarity, improving structural flow, and generating TikZ code for scientific diagrams. All diagrams are based on data collected and processed independently by the authors; the AI tools had no access to the underlying research data, ensuring the technical accuracy of the visualizations. Additionally, generative AI was used to produce the text-to-speech audio for the supplemental video. The authors reviewed and edited all outputs and take full responsibility for the content and integrity of the final publication and all associated multimedia. Acknowledgments This work was supported by the European Union under the ROBOPROX (Robotics and advanced industrial production) project (reg. no. CZ.02.01.01/00/22_008/0004590) and the Grant Agency of the Czech Technical University in Prague, grant No. SGS26/117/OHK3/1T/37. References Abreu, L. R. d., & Nagano, M. S. (2024). Contributions to new variants of the open shop scheduling problem: modeling and solution methods. Ph.D. thesis Universidade Federal do Ceará. Publisher: Galoá. Akbarzadeh, B., & Maenhout, B. (2024). A study on policy decisions to embed flexibility for reactive recovery in the planning and scheduling process in operating rooms. Omega, 126, 103061. doi:10.1016/ j.omega.2024.103061. Arkin, E. M., & Hassin, R. (1994). Approximation algorithms for the geometric covering salesman problem. Discrete Applied Mathematics, 55, 197–218. doi:10.1016/0166-218X(94)90008-6. Blank, J., & Deb, K. (2020). Pymoo: Multi-objective optimization in python. IEEE Access, 8, 89497–89509. doi:10.1109/access.2020.2990567. Carvalho, M., & Almeida, I. F. (2022). The role of pharmaceutical compounding in promoting medication adherence. Pharmaceuticals, 15, 1091. doi:10.3390/ph15091091. Dauzère-Pérès, S., Ding, J., Shen, L., & Tamssaouet, K. (2024). The flexible job shop scheduling problem: A review. European Journal of Operational Research, 314, 409–432. doi:https://doi.org/10.1016/j. ejor.2023.05.017. Deckerová, J., Váňa, P., & Faigl, J. (2024). Combinatorial lower bounds for the generalized traveling salesman problem with neighborhoods. Expert Systems with Applications, 258, 125185. doi:10.1016/j. eswa.2024.125185. 30 Fioravantes, F., Knop, D., Křištan, J. M., Melissinos, N., & Opler, M. (2024). Exact algorithms and lower- bounds for multiagent path finding: Power of treelike topology. In Proceedings of the AAAI Conference on Artificial Intelligence (p. 17380–17388). volume 38. Garey, M. R., & Johnson, D. S. (1979). Computers and intractability volume 174. freeman San Francisco. Go, M.-j. (2021). Lip factory by color tailor: Smart factory system. CES 2021 Innovation Awards Honoree, Amorepacific R&D Center. Hanzálek, Z., & Šucha, P. (2017). Time symmetry of resource constrained project scheduling with general temporal constraints and take-give resources. Annals of Operations Research, 248, 209–237. doi:10.1007/ s10479-016-2184-6. Hao, J., Rodriguez-Monguio, R., & Seoane-Vazquez, E. (2015). Fixed-dose combination drug approvals, patents and market exclusivities compared to single active ingredient pharmaceuticals. PLOS ONE, 10, e0140708. doi:10.1371/journal.pone.0140708. Hartmann, S., & Briskorn, D. (2022). An updated survey of variants and extensions of the resource- constrained project scheduling problem. European Journal of Operational Research, 297, 1–14. doi:10. 1016/j.ejor.2021.05.004. Heinz, V., Novák, A., Vlk, M., & Hanzálek, Z. (2022). Constraint programming and constructive heuristics for parallel machine scheduling with sequence-dependent setups and common servers. Computers & Industrial Engineering, 172, 108586. doi:10.1016/j.cie.2022.108586. Itai, A., Papadimitriou, C. H., & Szwarcfiter, J. L. (1982). Hamilton paths in grid graphs. SIAM Journal on Computing, 11, 676–686. doi:10.1137/0211056. Janning, K., Housin, A., Schulte, C., Erkens, F., Frenken, L., Herbst, L., Nießing, B., & Schmitt, R. (2025). Conflict-based model predictive control for multi-agent path finding-experimentally validated on a magnetic planar drive system. Frontiers in Control Engineering, 6, 1645918. doi:10.3389/fcteg.2025. 1645918. Laborie, P. (2003). Algorithms for propagating resource constraints in AI planning and scheduling: Existing approaches and new results. Artificial Intelligence, 143, 151–188. doi:10.1016/S0004-3702(02)00362-4. Leisch, F., Weingessel, A., & Hornik, K. (1998). On the generation of correlated artificial binary data. Work- ing Papers SFB "Adaptive Information Systems and Modelling in Economics and Management Science", . doi:10.57938/6884f809-93bc-4497-ab2b-fc1611198f5b. Maenhout, B., & Vanhoucke, M. (2013). An integrated nurse staffing and scheduling analysis for longer-term nursing staff allocation problems. Omega, 41, 485–499. doi:10.1016/j.omega.2012.01.002. National Center for Health Statistics (2012). National Health and Nutrition Examination Survey (NHANES) 2009–2010 Data Documentation, Codebook, and Frequencies: Prescription Medications (RXQ_RX_F). Centers for Disease Control and Prevention (CDC). URL: https://wwwn.cdc.gov/nchs/nhanes/ 2009-2010/RXQ_RX_F.htm accessed: 2026-01-22. 31 Nilsson, H., & Ternerot, J. (2022). Path planning algorithm for levitating planar motion system. Master’s thesis Lund University. Noon, C. E., & Bean, J. C. (1993). An efficient transformation of the generalized traveling salesman problem. INFOR: Information Systems and Operational Research, 31, 39–44. doi:10.1080/03155986. 1993.11732212. Quinn, K. J., & Shah, N. H. (2017). A dataset quantifying polypharmacy in the united states. Scientific Data, 4, 170167. doi:10.1038/sdata.2017.167. Rehme, L. (2025). Path-Planning Algorithms for a Levitating Planar Motion System. Master’s thesis Lund University. Sonntag, E., Stepanek, F., Vrba, J., Mutylo, E., & Slezakova, S. (2023). Method of manufacturing multiple dosage forms with personalised composition and device for carrying out this method. International Patent Application WO 2023/151732 A1. Stern, R. (2019). Multi-agent path finding – an overview. In G. S. Osipov, A. I. Panov, & K. S. Yakovlev (Eds.), Artificial Intelligence: 5th RAAI Summer School, Dolgoprudny, Russia, July 4–7, 2019, Tutorial Lectures (p. 96–115). Cham: Springer International Publishing. doi:10.1007/978-3-030-33274-7_6. Sutherland, J. J., Daly, T. M., Liu, X., Goldstein, K., Johnston, J. A., & Ryan, T. P. (2015). Co-prescription trends in a large cohort of subjects predict substantial drug-drug interactions. PLOS ONE, 10. doi:10. 1371/journal.pone.0118991. Vanhoorne, V., & Vervaet, C. (2020). Recent progress in continuous manufacturing of oral solid dosage forms. International Journal of Pharmaceutics, 579, 119194. doi:10.1016/j.ijpharm.2020.119194. Wang, Y., Ma, H.-S., Yang, J.-H., & Wang, K.-S. (2017). Industry 4.0: a way from mass customization to mass personalization production. Advances in Manufacturing, 5. doi:10.1007/s40436-017-0204-7. Wastesson, J. W., Morin, L., Tan, E. C. K., & Johnell, K. (2018). An update on the clinical consequences of polypharmacy in older adults: a narrative review. Expert Opinion on Drug Safety, 17, 1185–1196. doi:10.1080/14740338.2018.1546841. Yadav, A., & Jayswal, S. (2018). Modelling of flexible manufacturing system: a review. In- ternational Journal of Production Research, 56, 2464–2487. doi:10.1080/00207543.2017.1387302. arXiv:https://doi.org/10.1080/00207543.2017.1387302. Zahradka, D., Andreychuk, A., Kulich, M., & Yakovlev, K. (2022). Quality analysis of multi-agent multi- item pickup and delivery solutions using a decoupled approach. IFAC-PapersOnLine, 55, 61–66. doi:10. 1016/j.ifacol.2023.01.134. 32 Appendix A. Notation Symbol DescriptionSymbolDescription HTopology layout of the environment.n tiles Number of tiles. ISet of interface tiles.n inter Number of interfaces. n disp Number of dispensers.PSet of orders. MSet of movers.n movers Number of movers. m max Maximum number movers allowed.GSet of drugs. DPacked dispensersφ :D∪I →HPlacement function for dispensers to tiles. P history Set of historical order data.TSet of current orders. L tiles Distance matrix over the set of tiles.∆ g p Duration of dispensing of drug g for order p. S p Set of drugs requested by order p∈P.ASet of alternative operations. η ′ p Time to traverse the path for order p∈P.RSet of resting sites. η I Constant time for cartridge swap.WSet of transits between sites with idle periods. Table Appendix A.1: List of symbols. Appendix B. Packing: Maximization of Pairwise Correlations This section describes the additional step in the packing, which aims to maximize the sum of the pairwise correlations. The motivation comes from the correlation matrix in Figure Appendix B.1 that we derived from the real order data (National Center for Health Statistics, 2012). The matrix displays the pairwise correlation of the individual drugs appearing in the patients’ prescriptions. By manual investigation, we found that, for example, drugs for high blood pressure tend to co-occur with drugs for diabetes, i.e., a positive correlation. On the other hand, a negative correlation was discovered, for example, for a pair of drugs, both used for cardiac patients that work as alternatives to each other—one is a stronger variant of the other. 010203040 0 10 20 30 40 drug g [-] drug g [-] −0.2 0 0.2 Figure Appendix B.1: Example of drug correlation matrixO = o g,g ′ g,g ′ ∈G derived from real data (with subtracted main diagonal). Therefore, our idea would be to find, among all solutions with the minimal maximum tile utilization, the one that also maximizes the sum of pairwise correlations between drugs when placed on the same tile. In this way, we can reduce the number of steps, as often co-occurring drugs are filled on the same tile without 33 moving or transitioning to a different tile. Initially, we experimented with a modification of model (4.1)– (4.11) that would express this secondary objective via the scalarization technique. However, the results were largely unsatisfactory. Thus, we decided to split the solution into two phases: first, we solve model (4.1)–(4.11), and then fix the maximum tile utilization ˆμ max and the number of dispensers ˆz g for each drug g. After that, we solve a correlation model that aims to maximize pairwise correlations given the maximum utilization bound for each tile. The inputs to the correlation model are feasible values from the packing model (4.1)–(4.11), where ˆμ max is achieved maximum tile utilization, the number of tiles that drug g occupies ˆz g and the expected utilization of a single dispenser ˆπ g of drug g. These values are fixed in the correlation model and not subject to optimization; the only free variables are tile-drug assignments y k,g ∈0, 1: max n tiles X k=1 X g,g ′ ∈G,g̸=g ′ o g,g ′ · y k,g · y k,g ′ (B.1) subject to n tiles X k=1 y k,g = ˆz g ∀g ∈G(B.2) 1≤ X g∈G y k,g ≤ d max ∀k ∈1,...,n tiles (B.3) X g∈G ˆπ g · y k,g ≤ ˆμ max ∀k ∈1,...,n tiles − n inter (B.4) where y k,g ∈0, 1 ∀k ∈1,...,n tiles ,∀g ∈G.(B.5) The model (B.1)–(B.5) is an MIQP model, which resembles a variant of the quadratic multiple knapsack problem. The optimization of pairwise correlations indeed has some merits. For example, an example packing produced by the first stage has a total sum of pairwise correlations of −0.31, while after optimization, it is 1.88—with the same maximum tile utilization. Appendix C. Experiments: Extended Data Tables Some of the tables in the manuscript were condensed to highlight key trends and significant values for clarity. This section provides the extended versions of those tables. 34 Mover count Square C max Doubleline C max Ring C max Line C max SchedulingRoutingSchedulingRoutingSchedulingRoutingSchedulingRouting 133427+0 (0.0%)34506+0 (0.0%)35056+0 (0.0%)35233+0 (0.0%) 2 16722+7 (0.0%)17295+22 (0.1%)17556+14 (0.1%)17663+28 (0.2%) 311163+11 (0.1%)11524+21 (0.2%)11782+35 (0.3%)11788+31 (0.3%) 48373+10 (0.1%)8638+21 (0.2%)8806+29 (0.3%)8875+38 (0.4%) 5 6700+14 (0.2%)6943+23 (0.3%)7091+30 (0.4%)7131+39 (0.5%) 65602+17 (0.3%)5822+27 (0.5%)5947+41 (0.7%)5992+54 (0.9%) 74808+17 (0.4%)5030+24 (0.5%)5092+48 (0.9%)5151+75 (1.5%) 84224+24 (0.6%)4399+37 (0.8%)4512+46 (1.0%)4560+53 (1.2%) 93762+17 (0.5%)3961+33 (0.8%)4012+58 (1.4%)4137+72 (1.7%) 103388+19 (0.6%)3557+31 (0.9%)3623+73 (2.0%)3757+71 (1.9%) 11 3137+33 (1.1%)3310+50 (1.5%)3372+80 (2.4%)3659+104 (2.8%) 122942+23 (0.8%)3042+42 (1.4%)3108+56 (1.8%)3286+100 (3.0%) Table Appendix C.1: Full version of Table 4. Mover count Square C max Doubleline C max Ring C max Line C max SchedulingRoutingSchedulingRoutingSchedulingRoutingSchedulingRouting 14731+0 (0.0%)5731+0 (0.0%)6202+0 (0.0%)6201+0 (0.0%) 22386+5 (0.2%)2897+10 (0.3%)3181+13 (0.4%)3296+14 (0.4%) 31592+9 (0.6%)1980+21 (1.1%)2171+25 (1.2%)2167+32 (1.5%) 41199+10 (0.8%)1492+23 (1.5%)1714+42 (2.5%)1638+50 (3.1%) 5 962+16 (1.7%)1187+33 (2.8%)1362+27 (2.0%)1410+32 (2.3%) 6817+22 (2.7%)1012+43 (4.2%)1131+42 (3.7%)1158+45 (3.9%) 7710+22 (3.1%)935+36 (3.9%)969+46 (4.7%)987+44 (4.5%) 8619+23 (3.7%)858+39 (4.5%)914+52 (5.7%)980+53 (5.4%) 9 571+23 (4.0%)728+48 (6.6%)813+46 (5.7%)885+66 (7.5%) 10 524+24 (4.6%)720+50 (6.9%)743+41 (5.5%)834+69 (8.3%) 11500+34 (6.8%)642+58 (9.0%)712+43 (6.0%)807+107 (13.3%) 12 448+37 (8.3%)605+56 (9.3%)665+60 (9.0%)689+79 (11.5%) Table Appendix C.2: Full version of Table 5. 35