Energy-efficient Automated Vertical Farms

Autonomous vertical farms (VFs) are becoming increasingly more popular because they allow to grow food minimising water consumption and the use of pesticides, while greatly increasing the yield per square metre compared with traditional agriculture. To meet sustainability goals, however, VFs must operate at maximum efficiency; it would be otherwise impossible to compete with the energy source powering plant growth in traditional agriculture: the sun. We introduce the Vertical Farming Elevator Energy Minimisation Problem ( VFEEMP ), which arises when minimising the energy consumption of automatic elevators servicing VFs. We prove that the decision problem associated with the VFEEMP is NP -complete. To solve the problem, we propose three Mixed-Integer Linear Programming (MIP) formulations together with valid inequalities, and a Constraint Programming model. We present a large set of instances, both synthetic and derived from real-life data, and we determine through extensive computational experiments which instance characteristics have an impact on the difficulty of the problem and which formulations are the most suitable to solve the VFEEMP . demand with the available infrastructure), the number of tray movements (trays are allowed to move from a shelf to another one during their growth), the number of shelf reconfigurations (i.e., the number of times the operator changes growth conditions for a shelf), and the number of shelves used. They prove that the problem is NP -complete and provide MIP models for the four versions. The models are tested on instances with up to 6 crops, 12 shelves, and a time horizon of 100 days. Their results show that the choice of the objective function heavily influences both the computational performance of the models and the characteristics of the solution, hinting that a multi-objective approach could be useful to build a real-world decision support system.


Introduction
Access to high-quality food produced sustainably is a pressing problem in both developed and developing countries. Standard agricultural practices are not well-suited to scale for a world with a larger and increasingly urbanised population, less available water, and climate change [8]. Simply increasing the amount of land dedicated to agriculture is not a compelling long-term option because of the associated side-effects: deforestation, soil depletion, the need to use pesticides and fertilisers, emissions due to transport between the place of production and that of consumption, etc. Moreover, standard agriculture provides low yields in terms of the amount of nutrients produced per square meter of land used [4]. To mitigate the effects that the forecast increase of food demand will have on urban and rural communities, significant efforts are devoted to devise alternative agricultural practices. Among the main objectives pursued is to increase the yield per square meter and reduce externalities on the environment, while staying economically competitive [3].
Vertical Farming (VF), an agricultural technique which is gaining increasing traction, is establishing itself as an invaluable tool to meet the above objectives. VF consists of growing crops in stacked layers hosted in indoor support structures. The hosting infrastructure provides the plants with everything they need to grow in optimal conditions: the right amount of light and water, ventilation and appropriate CO 2 levels, nutrients, protection from pests, controlled temperatures. Such infrastructures vary considerably in terms of complexity, level of control they offer over the growth environment, and size [13,12].
All types of structures, irrespective of their size or level of sophistication, provide some basic advantages compared to traditional farming: they protect the plants from weather and climate variability; they shelter them from parasites and pathogen; they allow to control the water and soil to ensure it is not affected by the presence of heavy metals or other dangerous substances; they can be employed everywhere, thus bringing production closer to consumption and reducing transportation costs and emissions. For example, supermarkets, restaurants and hotels are increasingly adopting VF cabinets, approximately as large as standard refrigerators, to grow in-house herbs, leafy vegetables and berries. However, for VF to make a considerable impact towards achieving a more sustainable food supply chain and to be economically viable, operators strive to operate at larger scales [4]. Structures of the size of silos or re-purposed city buildings are more appropriate to host crops central to human diets, such as staple foods, vegetables and fruits [7,2]. Although the technology which would allow such very large-scale projects to operate is not yet fully developed, interesting concepts are starting to emerge. For example, 6-to-12-metre tall "growth towers" in which plants grow under completely controlled conditions and with minimal human interaction [9] are available on the market. In such towers, each tray hosting a crop receives light, water, nutrients and ventilation from a computer-controlled Internet-of-Things system. Trays are moved around, irrigated, and inspected via automatic elevators equipped with water hoses and cameras. A human operator is needed to bring a new tray to the elevator's loading bay, and to collect it at the end of its growth period.
These large-scale systems allow a lower energetic footprint. For example, they are large enough to be equipped with solar panels for electricity production and contain enough air that water can be harvested from condensed humidity. In other respects, however, VF is still more energy-consuming than traditional farming. First, it uses artificial light to provide crops with an optimal amount of lighting and faster growth times. Although new technology such as LED lighting consumes significantly less energy compared to systems from just 5 years ago, it is clearly not as efficient as using sunlight. Second, for large-scale systems such as the growth towers described above, the automated elevators which serve the tower use a significant portion of the energy required by the structure.
The objective of this paper is to develop mathematical models and decision support tools to lower the energy consumption of the automatic elevators servicing large-scale VF structures. To do so, in the rest of this section, we give a precise description of the problem. Section 3 presents existing literature on the optimisation of vertical farms and on other related problems. We formalise our problem and prove its N P-completeness in Section 4. In Section 5 we present three MIP formulations and a number of valid inequalities aimed at strengthening their linear relaxations, and we also present a formulation based on Constraint Programming (CP). After describing in Section 6 the instance sets we use, we present the results of a large set of computational experiments in Section 7. We conclude and point out future research directions in Section 8.
The main contributions of this paper are the following: • We introduce a problem arising from the operation of autonomous vertical farms, which has timely applications in a rapidly growing industry. • We prove that the problem is N P-complete by reduction from the well-known SubsetSum problem. • We present three MIP formulations and strengthen them with valid inequalities and variable fixing techniques, as well as a CP formulation. • We generate and publish a large set of instances, both synthetic and derived from real-life data. We also determine which instance characteristics have an impact on the difficulty of the problem. • We identify the best-performing formulation and set of valid inequalities. We also make public the source code of all our solvers.

Problem description
A VF tower is composed of stacked shelves in which a planner places trays containing crops. Each tray will spend some time in this system until the crops it contains are ready to be harvested. Each shelf can host one tray at a time and, after a tray leaves a shelf, another one can take its place. Trays become available from a "depot" at the bottom, spend the required amount of time on a shelf and, at the end of the crop's growth period, return to the depot. In an automatic vertical farm, an elevator moves the trays between the depot and the shelves. Figure 1 gives a schematic representation of a stack with five shelves (plus the depot) served by an elevator. A human operator brings a seeded tray to the depot so that the elevator can place it into a shelf. At the end of the growth cycle of the crop, the elevator brings the tray down and the operator picks it up from the depot and moves it into an area where the crop will be harvested. The elevator also serves the trays performing activities such as watering and monitoring. For example, it can take pictures that a computer vision algorithm processes to ensure crops are growing well and are not affected by pests. Each tray has an associated list of tasks, and each task requires that the elevator travels to the corresponding shelf and spends some time there.
Each tray has a start time window during which it must be picked up at the depot to be placed on its assigned shelf. We define the time at which the elevator picks up the tray to bring it to its shelf as the tray's start time. There are also time windows related to each task and to the final pick up of the tray, but these are relative to the tray's start time. For example, if a tray has a start time window of one day (from hour 0 to 24), it can be brought to its shelf at hour 12. Then, if it must be watered after 6-8 hours from the growth start, this would translate into a task time window spanning from hour 18 to hour 20.
A computer system is in charge of deciding in which order the tasks should be completed, and correspondingly moving the elevator to perform each task. The objective of the system is to minimise the energy consumption of the elevator, making sure that all tasks are performed within their time windows. The elevator uses energy to travel between shelves. Therefore, the energy minimisation objective translates into the minimisation of the distance travelled by the elevator.
In the problem considered in this paper, another subsystem in the control software of the tower has already taken care of assigning each tray to the shelf it will occupy (see Section 8 for possible extensions). The elevator must then transport each tray to and from its assigned shelf, and visit the trays to perform the given list of tasks. The elevator can carry out tasks associated with different trays in any order, as long as the time windows are respected. On the other hand, tasks associated with the same tray must be carried out in the given order.
The objective of the problem is to determine a sequence of operations which minimises the travel distance of the elevator. Because the travel time of the elevator (a few seconds) is negligible compared to the time it must spend at the shelves to perform the tasks (several minutes), we can disregard it when timing the operations.

Literature review
There is a limited but growing amount of literature on the topic of the optimisation of VF. Two of the earliest works, by Yang, Hari, and Kuo [20] and Bennell, Martinez, and Potts [5], respectively, are only available as short abstracts in conference proceedings.
Yang, Hari, and Kuo [20] study the problem of selecting and scheduling a mix of crops to grow in a VF system at a tactical level. The objective is to maximise profits, taking into account that the price of each crop fluctuates during the year. The fluctuation is due to the possibility of growing some of the crops using open-air agriculture, which is highly seasonal. Therefore, it is more profitable to harvest a crop in the VF system during the period of the year when it is not available using open-air agriculture. The author proposes a small case study with a 2-year time horizon and 13 possible crops, which is solved using a MIP.
The problem of meeting as much demand as possible using crops grown in a VF cabinet was studied by Bennell, Martinez, and Potts [5]. The authors consider a stack of shelves similar to the one presented in Figure 1, although smaller in scale: their proposed use-case is a cabinet placed directly at the consumer site (e.g., at a restaurant or a supermarket). In their settings, the capacity of each shelf and the growth speed of each crop are determined by a set of growth conditions, which the operator can change. For example, each crop can grow on a different growth medium, receive more or less light, water, etc. Using these factors to influence the length of the growth cycles, Bennell, Martinez, and Potts minimise the amount of unmet demand, weighted by crop prices. The authors propose a heuristic rolling-horizon algorithm.
Yang, Huang, and Ang [21] expand on [20], considering again the problem of selecting and scheduling crop growth in a VF system. The authors consider a setting in which crops receive sunlight (as opposed to artificial light) and the amount of light depends on the shelf's position and on whether the system operator installs rooftop solar panels, which obstruct the glass ceiling. Similar to [20], the authors aim at maximising the profit considering the seasonal variability of crop prices. Additionally, they introduce crop adjacency considerations, which make it more desirable for crops of the same family to grow on adjacent shelves. The MIP formulation proposed by the authors is not able to scale to realistically sized instances and, therefore, they implement a greedy heuristic which provides solutions for time horizons of up to 50 months.
Recently, Santini et al. [17] considered the problem of planning the production of crops on VF cabinets. The authors consider a system in which crops grow under controlled conditions (soil, water, light, CO 2 levels, etc.), which can change day-by-day and on a shelf-by-shelf level. They propose four objective functions which a planner might want to minimise: the amount of unmet demand (if it is not possible to satisfy the entire demand with the available infrastructure), the number of tray movements (trays are allowed to move from a shelf to another one during their growth), the number of shelf reconfigurations (i.e., the number of times the operator changes growth conditions for a shelf), and the number of shelves used. They prove that the problem is N P-complete and provide MIP models for the four versions. The models are tested on instances with up to 6 crops, 12 shelves, and a time horizon of 100 days. Their results show that the choice of the objective function heavily influences both the computational performance of the models and the characteristics of the solution, hinting that a multi-objective approach could be useful to build a real-world decision support system.
In another recent paper, Cetegen and Stuber [6] present a robust optimisation model to design a VF system and schedule crop growth. The authors aim at building a portfolio of crops to grow, to minimise the risk of economic loss. Because controlled-environment agricultural techniques such as VF virtually eliminate risks on the grower's side (bad weather, pests, etc.) the authors consider risks coming from the market's side, i.e., demand and price variability. Using a cutting-plane method to solve a semi-infinite programme with semi-infinite constraints, the authors apply their algorithm to two case studies, validating the economic viability of using VF to grow crops such as lettuce, spinach, tomatoes, strawberries and mushrooms.
We also highlight a recent Master's thesis [1], which deals with task scheduling in a semiautomatic VF plant factory. The objective considered is the minimisation of the makespan of tasks to be performed by sensors and actuators on the crop trays growing in the VF system. This objective allows the authors to model the problem as a Flexible Job-Shop Scheduling Problem and use a variant of the Genetic Algorithm by Pezzella, Morganti, and Ciaschetti [14] to solve it.
In the broader context of reducing the environmental footprint of industrial production processes, we place our work within a growing number of papers which explicitly consider sustainability and energy efficiency when optimising production planning. For example, Wu and Che [19] propose a bi-objective flow-shop scheduling problem in which the two objectives are the minimisation of the total make-span, and the minimisation of energy consumption. Qiu, Qiao, and Pardalos [15] present a production routing problem in which excessive carbon emissions are penalised in the objective function, while reduced emissions yield a prize.
Finally, we remark that our problem has a superficial similarity with the Travelling Salesman Problem with Time Windows (TSPTW): if the task time windows were absolute rather than relative to the planting time, our problem could be seen as an extension of a 1-dimensional TSPTW. The current literature, however, does not present any work neither on routing problems in which the time window of a customer is relative to the visit time of another customer, nor on 1-dimensional special cases of routing problems with time windows.

Problem formalisation
Consider a tower with n shelves, making up set S = {1, . . . , n}. We denote the bottom depot as 0 and the set of all shelves plus the depot as S ′ = {0} ∪ S. During the planning horizon, m trays will be placed on the shelves; we denote them as T = {1, . . . , m}. The time horizon itself is discretised in time intervals and denoted as Ω = {1, . . . , ω}.
Each tray t ∈ T has a start time window (withw 0 t ,w 0 t ∈ Ω andw 0 t ≤w 0 t ) during which it must be picked up at the depot. The elevator must stay at the depot for d 0 t ≥ 1 additional time units when picking up the tray, to allow loading operations. The start time window is absolute. Conversely, the other time windows introduced below will be relative to the time in which the elevator reached the depot to pick up the tray. Furthermore, there is a list of l t tasks which the elevator must perform on tray t; the set of these tasks is denoted as J t = {1, . . . , l t }. Each task j ∈ J t has an associated (relative) time window during which the elevator must start the task, and a duration d j t ≥ 1 during which the elevator cannot move away from the shelf containing t. Finally, we denote with the (relative) time window during which the elevator must pick up tray t and bring it back to the depot, and with d l1+1 t ≥ 1 the additional time the elevator must spend at the depot after bringing the tray, to allow unloading operations. We define the set of extended tasks as . We allow time windows to overlap for a given tray, i.e., it could happen thatw j t +d j t ≥w j+1 t , ∀j ∈ {1, . . . , l t }. However, to ensure the proper growth of the crops, we require tasks to be executed in their given order, i.e., for a given tray, a task j + 1 cannot be performed before task j.
It is also possible that the allowed "shelf life" of a tray (i.e., the moments in which it is allowed to be present on its assigned shelf) can overlap with that of another tray assigned to the same shelf. In other words, given two trays t 1 , t 2 such that both are assigned to the same shelf, it can happen thatw 0 t1 +w However, the two intervals [w 0 t1 ,w 0 t1 +w ] (if non-empty) cannot overlap, as they represent the time instants in which the trays are required to be on the shelf and two trays cannot occupy the same shelf at the same time.
Because all shelves are equally spaced in the silo, the distance travelled by the elevator when moving from a shelf s 1 ∈ S ′ to another shelf s 2 ∈ S ′ can be easily expressed as c s1,s2 = c s2,s1 = |s 1 − s 2 |. We also denote the shelf assigned to a tray t ∈ T as s(t) ∈ S.
We define the Vertical Farming Elevator Energy Minimisation Problem (VFEEMP) as the problem of deciding in which order to perform the tasks, to minimise the distance travelled by the elevator. Differently from any other problem we are aware of, most time windows are relative and depend on other decisions taken by the planner, namely the start time of the planting task. Because of this distinguishing feature, the VFEEMP seems to be neither a special case nor a generalisation of other known optimisation problems and, thus, a first interesting question regards its computational complexity. Theorem 1 answers this question, proving that the decision version of the VFEEMP is N P-complete.

Theorem 1. Let VFEEMP-Decision be the decision problem which asks whether there is a solution to the VFEEMP of cost at most x, for a given
Proof. We prove the theorem by proving a stronger statement: that determining whether VFEEMP is feasible is itself N P-complete. This clearly implies that VFEEMP-Decision is N P-complete.
Consider an instance of the N P-complete SubsetSum problem: given a set of n positive integers B = {b 1 , . . . , b n } and a positive integer c, the problem asks to determine whether there is a subset S ⊆ B such that ∑ b∈S b = c. We show how to create an instance of the VFEEMP such that, if the instance is feasible, then SubsetSum has answer Yes and, otherwise, it has answer No.
The main idea is to first transform the SubsetSum instance by multiplying each b i and c by 2, to only deal with even numbers. Next, we build an instance in which each tray uses 2b i units of elevator time, in which the elevator can never be idle, and in which there is a special tray which must be picked up at time 2c + 1. In this way, being able to serve a subset of trays before time 2c + 1 corresponds to finding a subset of the (doubled) integers summing up to 2c, i.e., a set of the original integers which sums to c.
Formally, we create the VFEEMP instance as follows: let B = 2 ∑ n i=1 b i , set the number of shelves to 1, the number of trays to n + 1 and the time horizon to Ω = {1, . . . , B + 4}. Each of the first i = 1, . . . , n trays has: • Start time window  Figure 2: Visualisation of a feasible VFEEMP transformed from SubsetSum.
meaning that the end task must be processed right after the start task.. Tray n + 1 has: • Start time window . In our instance this relative time window is equivalent to absolute time window and d 2 n+1 = 1. In our instance this relative time window is equivalent to absolute time window If VFEEMP is feasible, let (π 1 , . . . , π n+1 ) be the permutation corresponding to a solution and let j be the index such that π j = n + 1. Then, because each tray i consumes exactly 2b i units of elevator time and because tray n + 1 must be picked up exactly at time 2c + 1, then we must have that SubsetSum has the answer Yes. Figure 2 depicts such a situation.
Noting that the above VFEEMP instance can be built in linear time in the size n of the SubsetSum instance concludes the proof.

Model formulation
In this section, we propose four formulations for the VFEEMP. The first three are based on MIP: the first uses a polynomial number of variables and disjunctive constraints; the second one uses a pseudo-polynomial number of variables; the third one is a hybrid which only uses binary variables from the first and the second formulation. Finally, the fourth formulation is based on CP.

A model with disjunctive constraints
To model the fact that the elevator starts and must return to the depot, we first add to Ω a dummy time interval 0, and to T a dummy tray with index 0, assigned to the depot with start time window [ 0, 0 ] and d 0 0 = 1, no intermediary tasks, and end time window . Variable y t1,j1,t2,j2 takes value 1 iff the elevator performs task j 1 on tray t 1 and, immediately afterwards, it performs task j 2 on tray t 2 . Recall that the assignment of trays to shelves is given, and that we denote with s(t) ∈ S the shelf assigned to tray t ∈ T . We denote with c(y t1,j1,t2,j2 ) ∈ N the travel cost associated with variable y t1,j1,t2,j2 . To assign a value to c(·), we consider the following cases: Figure 3: Possible cases for the costs associated with variables y t 1 ,j 1 ,t 2 ,j 2 .
The elevator movements associated with each of these cases are shown in Figure 3. Note that a case in the figure depicts more movements than those considered in the corresponding cost to prevent counting some movement costs twice. For example, in Case (10) we do not count the cost of bringing back the elevator to the depot after picking up tray t 1 because it is already counted in the previous move, which was either a Case (4), a Case (7), or another Case (10). We also introduce the following notation: In other words, α tj and β tj are the earliest and latest time intervals at which task j ∈ J ′ t can start. We now notice that not all possible quadruples (t 1 , j 1 , t 2 , j 2 ) are valid indices for variables y: • If t 2 = 0 and j 2 = 0, then no indices t 1 , j 1 , t 2 , j 2 can be valid, as the start task of the dummy tray cannot have any predecessor. • Analogously, if t 1 = 0 and j 1 = 1, then no indices t 1 , j 1 , t 2 , j 2 can be valid, because the end task of the dummy tray cannot have any successor. • We can also omit indices which would result in time window violations. Specifically, if β t2j2 < α t1j1 + d j1 t1 then task j 2 must be performed before task j 1 and cannot be its successor. Therefore, we need not consider indices t 1 , j 1 , t 2 , j 2 satisfying the above condition.
• If t 1 = t 2 , then only indices for which j 2 = j 1 + 1 are valid; otherwise, tasks would be performed in the wrong order for that tray. • If there is another task j 3 ∈ J ′ t3 (for some tray t 3 ∈ T ) which must take place between j 1 and j 2 , then these two tasks cannot happen in succession. More precisely, indices t 1 , j 1 , t 2 , j 2 are not valid if Let us denote with V the set of valid indices for variables y, i.e., the set of all quadruples (t 1 , j 1 , t 2 , j 2 ) which are not excluded by one of the above conditions. We also consider continuous variables s tj indicating the starting time of task j on tray t (t ∈ T, j ∈ J ′ t ). Finally, let ∆ t1t2 ∈ {0, 1} be a binary variable defined for each pair of trays t 1 , t 2 ∈ T such that s(t 1 ) = s(t 2 ). It will take value 1 iff tray t 1 is placed on the shelf before tray t 2 . The VFEEMP can then be modelled with the following MIP, which we denote M1: The objective function (13) minimises the travel cost of the elevator. Constraints (14) ensure that every task of every tray (excluding the start task of the dummy tray) has a predecessor and, therefore, is processed. Constraints (15) force the elevator to process the start task of the dummy tray, and constraints (16) make sure that every task (excepted those associated with the dummy tray) has a direct predecessor and a direct successor. Constraints (17) and (18) make sure that each tray is picked up during its associated (absolute) start time window, and constraints (19) and (20) ensure that every task other than the start task takes place during its associated (relative) time window. Constraints (21) force tasks associated with the same tray to be scheduled in the given order and are necessary because task time windows can overlap. Constraints (22), in which M is a sufficiently large number, force the starting time of any task to be not earlier than the starting time plus the duration of its predecessor. The tightest value for M is β t1j1 + d j1 t1 − α t2j2 . Note that constraints (22) could easily be updated to incorporate a travelling time between shelves s(t 1 ) and s(t 2 ) if the time spent in elevator movements was not negligible compared to the time spent performing the tasks. Constraints similar to (22) are also used in the literature on routing problems with time windows to sequence customer visits (see, e.g., constraint (14) of [11]) and in scheduling to sequence tasks (see, e.g., constraint (5) of [18]). Constraints (23) link variables s and ∆ (again, using a sufficiently large number M ), and constraints (24) ensure that no two trays occupy the same shelf at the same time. The tightest value for M in (23) is β t2,lt 2 +1 + d lt 2 +1 t2 − α t10 . Variables ∆ and constraints (23)-(24) are required because tray shelf lives can overlap and, therefore, the model must determine a processing order for trays on the same shelf. Finally, (25), (26), and (27) are variable domain definitions.
For inequalities (22) and (23) we use big-M constraints over SOS/indicator constraints because, experimentally, the M factors do not cause numerical difficulties. Indeed, we use constraint-specific M values resulting in coefficients that are not much larger than the other coefficients in the model.

Valid inequalities for M1
We introduce the following sets of valid inequalities for model M1. Their impact is tested by means of computational experiments presented in Section 7.1.

Task incompatibility (TaskInc).
The task incompatibility considerations we use to reduce the set of feasible indices for variables y can be extended to more than two tasks at a time. For example, the following inequality is valid for model M1: In this case, (28) states that tasks j 1 , j 2 and j 3 cannot be performed in the given order, due to their time windows and durations. Although one could consider similar clique inequalities for larger sets of mutually incompatible tasks, detecting such sets becomes harder while the constraints would be less likely to be active. For this reason, we only add triangle inequalities of type (28) to model M1.
Force task sequence (TaskSeq). While constraint (28) excludes incompatible tasks, we next present a constraint aimed at forcing tasks to happen one immediately after the other. In other words, we look for suitable conditions on (t 1 , j 1 , t 2 , j 2 ) ∈ V which force y t1,j1,t2,j2 = 1. These conditions must ensure that (i) j 2 cannot start before j 1 , and that (ii) no other task j 3 (relative to some tray t 3 ) can start between j 1 and j 2 . A sufficient condition for (i) is that starting j 2 at its earliest possible time would make the elevator miss j 1 's time window. A sufficient condition for (ii) is that there is no valid start time for j 3 which is compatible with both j 1 's and j 2 's time windows. According to these considerations, we introduce the following valid variable fixing equality for model M1: Furthermore, if condition (ii) holds but condition (i) does not, then no task j 3 can start between tasks j 1 and j 2 ; however any of j 1 and j 2 can precede the other. In this case, if condition (ii) also holds when swapping j 1 and j 2 (that is, if no task j 3 can start between j 2 and j 1 , as well), then we can add the following constraint: Constraint (30) states that either j 1 is performed immediately before j 2 , or j 2 is performed immediately after j 1 . It thus excludes the possibility that the elevator performs the third task between j 1 and j 2 , no matter in which order j 1 and j 2 are performed.

2-and 3-cycle elimination (CycElim).
We next propose two families of valid inequalities similar to cycle elimination constraints in routing and elementary shortest path problems [see, e.g., 10]: Constraint (31), a 2-cycle elimination constraint, states that either j 1 takes place before j 2 , or j 2 takes place before j 1 . Constraint (32), a 3-cycle elimination constraint, extends the same reasoning to three tasks at a time. It forbids solutions inducing a "cycle" with j 1 performed before j 2 , j 2 before j 3 , and j 3 before j 1 .
Bounds on task start times (STBound). The next family of valid inequalities aims at bounding variables s tj both from above and below. Given a task j ∈ J ′ t , we first consider the set Π tj of tray-tasks pairs containing tasks which the elevator cannot perform after j and must, therefore, precede j: denotes the earliest possible finish time of task j ′ . Taking the maximum such value among tasks j ′ ∈ Π tj ensures that the elevator cannot be free to perform task j before this time. Therefore, the following inequality is valid: Analogously, we can consider the set of tray-task pairs which the elevator cannot perform before j and must, therefore, follow j: Then, the following inequality, analogous to (33), is valid: Minimum inter-task time (MinIT). If a task j ′ ∈ J t ′ must follow another task j ∈ J t (i.e., if (j ′ , t ′ ) ∈ Σ jt ) then, j ′ cannot start before s jt + d j t . Therefore, the following constraint is valid: Note that inequality (35) is not implied by (22) and (33). Furthermore, in the continuous relaxation of the problem, fractional values of y can produce a solution in which (34) is satisfied, but (35) is violated.

A model with a pseudo-polynomial number of variables
We introduce parameter s(t, j) for t ∈ T and j ∈ J ′ t , which represents the shelf where the elevator is located after performing task j: We now define a binary variable x tjk for t ∈ T , j ∈ J ′ t and α tj ≤ k ≤ β tj . Variable x tjk takes value 1 iff the elevator starts performing task j on tray t at time k. We also define an integer variable z k (k ∈ Ω \ {1}) indicating the cost of the movements made between time k − 1 and time k. Finally, we introduce a binary variable p ik (i ∈ S ′ and k ∈ Ω) taking value 1 iff the elevator is at shelf i at time k.
Then model M2 reads as follows: z k ≥ 0 and integer ∀k ∈ Ω, where Ω tjk contains the starting time indexes of task j of tray t for which the task would not be completed by time k. The objective function (36) minimises the travel cost of the elevator. Constraints (37) ensure that every task of every tray is processed (exactly once), and constraints (38) make sure that the elevator is at a unique position at any moment. Constraints (39) forbid the elevator to process more than one task at once. Constraints (40) force the elevator to be at the shelf assigned to the tray whose task it is starting to process. These constraints, together with (39) and objective function (36) ensure that the elevator will not move for the entire duration of the task. Constraints (41) ensure that each task (except the start task) is started during its relative time window. These constraints also implicitly ensure that the start task begins within its absolute time window. Constraints (42) force the tasks to be scheduled in the given order, and are required because the time windows of consecutive tasks can overlap. Constraints (43) make sure that the starting task of a tray can only take place once its associated shelf is free. We need these constraints because the allowed shelf lives can overlap for trays assigned to the same shelf. Finally, constraints (44) and (45)

Valid inequalities for M2
We adapt the families of valid inequalities STBound and MinIT devised for M1 (see Section 5.1.1). Valid inequalities (33) and (34) result in fixing x variables to 0 in M2: Inequality (35) corresponds to the following constraint:

A hybrid model
We now introduce a model which uses both binary variables y t1,j1,t2,j2 of M1 and x tjk of M2, and does not need any other set of variables. We denote this model as M3.
min ∑ (t1,j1,t2,j2)∈V c(y t1,j1,t2,j2 ) · y t1,j1,t2,j2 In the formulation of M3 we used the following sets: The objective function (52) and constraints (53)-(55) are the same as, respectively, objective (13) and constraints (14) (61) and (62), which link variables x and y. In particular, (61) prevents the elevator from starting a task j 1 at time k 1 if it is the predecessor of another task j 2 that has to be started at time k 1 + d j1 t1 − 1 or before. Analogously, (62) forbids j 1 to start at time k 1 if it is the predecessor of another task j 2 starting at a time k 2 such that k 1 + d j1 t1 > k 2 . Inequalities (28)-(32) and (49)-(51) are also valid for M3. Note that if the travelling time between shelves is not negligible, we only have to change (61) by adjusting the range of coefficient k 2 in the condition "k 2 < k 1 + d j1 t1 ".

A constraint programming model
Our CP model, M4, uses two sets of variables. The first set is comprised of integer variables s tj ∈ [ α tj , β tj ] indicating the starting time of task j ∈ J ′ t of tray t ∈ T . The second set is made of interval sequence variables A(t, j), with A(t 1 , j 1 ) = (t 2 , j 2 ) iff the elevator performs task j 2 on tray t 2 immediately after performing task j 1 on tray t 1 . The model reads as follows: The objective function (65) minimises the travel cost of the elevator by using the interval sequence variables to compute the appropriate values of function c(·) introduced in Section 5.1.
Interval sequence variables also ensure that all but one task (the last of the sequence) have a direct successor and that every task is processed. Constraints (66) and (67) make sure that every task (except the start tasks) takes place during its associated (relative) time window. Constraints (68) force tasks associated to the same tray to be scheduled in the given order. Constraint (69) forbids the elevator to process more than one task at a given moment. Constraints (70) ensure that no two trays occupy the same shelf at the same time. Domain constraints (71) and (72) define the variable domains, with (71) also ensuring that every starting task is processed during its absolute time window. Implementation details. We solve model M4 with the Cplex constraint programming optimiser (see Section 7). We define each task as an IloIntervalVar variable and model the domain constraints with functions setStartMin and setStartMax. Precedence constraints between the starting task and the other tasks of the same tray are enforced using IloStartBeforeStart and IloStartOf: the former to forbid that task start too early and the latter to forbid that they start too late. We model precedence constraints between two consecutive tasks of the same tray using IloEndBeforeStart. We used IloStartOf, IloEndOf, and the logical operator "||" to model constraints (70). We embed the IloIntervalVar variables in an IloIntervalSequenceVar variable, which we use to compute the objective function, with the help of the IloTypeOfNext function. Finally, we model the non-overlap constraint using IloNoOverlap.

Instance generation
We generated two sets of instances, named Synthetic and Realistic. The Synthetic set provides a wide variety of instance parameters (time horizon length, number of shelves and trays, number and duration of tasks, etc.) in order to test the quality of the proposed formulations on a diverse test body. The Realistic set provides instances which closely resemble real-life applications in VF silos.

Synthetic Instances
We describe the process used to generate feasible synthetic instances. Note that, because the VFEEMP is N P-complete, verifying the feasibility of an instance is, in general, hard. Therefore, we approach the instance generation problem starting from a backbone instance which is guaranteed to be feasible and transforming it into a complete instance while preserving feasibility. We use three phases to generate an instance: • We first generate a sequence of tasks with their respective start times and durations as they would appear in a feasible solution, i.e., we guarantee that the elevator can perform all tasks. We call this structure the backbone instance. • Next, we assign the tasks to the trays and the trays to shelves, ensuring that no two trays occupy the same shelf at the same time. • Finally, we add time windows around the start times generated to allow for optimisation opportunities and reordering of tasks.

Generating task start times and durations.
Let m ∈ N be the number of trays in the instance, ω ∈ N be the time horizon length, and δ ∈ N be the target number of tasks to assign to each tray. We first create a list of η := mδ tasks specifying, for each of them, their start time and their duration. Assigning tasks to trays. We first assign the start task of the trays; afterwards, we assign all other tasks. For each tray t ∈ [ 1, m ] , we select a task uniformly at random as t's start task. In this phase we exclude from the possible tasks those whose start time is too close to the end of the time horizon; to this end, we do not consider tasks j such that k s j > 4 5 ω. We remove all tasks assigned as start tasks from the list of available tasks. Next, we assign δ − 1 further tasks to each tray. Let j s t ∈ {1, . . . , η} be the index of the task chosen as start task for t and recall that tasks are indexed in increasing order of their start time. For each tray t, we scan available tasks one by one starting from j s t + 1 (or whatever is the first available task after j s t ). When scanning a task, we assign it to tray t with a probability π gr ∈ (0, 1). In case the task is assigned, we remove it from the list of assigned tasks; otherwise, we move on to scanning the next available task until we add all δ − 1 remaining tasks to tray t. Probability π gr measures how "greedily" we assign consecutive tasks to the same tray. A high value of π gr gives trays with shorter and non-overlapping shelf lives; conversely, low values of π gr result in trays with longer and overlapping shelf lives.
Assigning trays to shelves. We assign trays to shelves greedily. We first sort trays by their shelf life start and then assign each tray to the available shelf with the lowest index. A shelf is available if it is not occupied by another tray at that moment. The number of shelves is not bounded a priori, but we note that high values of π gr allow more trays to occupy the same shelf (during different periods, of course) and result in a number of shelves typically between m 2 and 3m 4 in our instances. Low values of π gr reduce the probability that two trays can be assigned to the same shelf and result in a number of shelves typically between 3m 4 and m.
Generating time windows. The last phase in our instance generation procedure consists of generating time windows around the task start times of the backbone instance. Given a task j, and time windows for non-start tasks are adjusted compared to the tray planting time to make them relative. The time window half-width is defined as ξ := γ ω mδ , where γ ≥ 1 is a parameter. Larger values of γ result in larger and more frequently overlapping time windows, which give more feasible task orderings for the elevator. Table 1 describes the parameters used for instance generation. We generate an instance for each possible combination of the parameters, which vary independently. This gives a total of 3 5 = 243 instances in the Synthetic set.  Figure 4 shows the optimal solution of an example instance generated with parameters m = 16, ω = 250, δ = 5, π gr = 0.6, γ = 1.5. The x axis represents the time instants and the y axis the shelf number. Each yellow box denotes the "shelf life" of a tray, whose progressive id number is reported next to the top-left corner of the box. The lines represent the elevator moving and performing tasks (in black) or being idle (in red).

Realistic instances
To generate the Realistic set, we used real-life confidential data about the growth cycle of six crops and the relative tasks to be performed when they grow in a VF silo. Compared to the Synthetic set, these instances have much larger time horizons, larger shelf lives and time windows, and more tasks per tray. Time horizons range from 15 to 128 days which, at a resolution of one time instant per 15 minutes, corresponds to values of ω between 1440 and 12 288. The instances contain either 5, 10, 15, or 20 shelves and between 5 and 40 trays. Each tray has, on average, 115.4 assigned tasks (including the start and end tasks).
The Realistic instances are particularly challenging because the number of variables and constraints in our models grows as a function of either the number of trays, the number of tasks, the time horizon length, the time windows sizes, or a combination of these factors. Furthermore, large time windows with significant overlaps contribute to the difficulty of these instances (see Section 7.1.3). Figure 5 shows the optimal solution of a Realistic instance with six trays and a time horizon of 80 days. Comparing Figure 4 and Figure 5, we highlight the differences in the length of the time horizon, the average shelf life, and the number of tasks associated with each tray.

Computational results
In this section, we present insights on the behaviour of our models on both the Synthetic and Realistic instances presented in Section 6. We analyse the impact of the valid inequalities introduced in Section 5 and identify instance characteristics which make the problem harder to solve. All models were implemented in C++; we solved the MIP models using Gurobi version 9 and the CP model using Cplex version 12.1. We make available on GitHub [16] the instance generator, the solver source code and the instance files. We performed the computational experiments on a 4-core machine equipped with an Intel Xeon CPU runing at 2.4GHz and 4GB RAM.

Results on the Synthetic instances
We used the synthetic instances as a test bed to investigate: • The impact of valid inequalities TaskInc, TaskSeq, CycElim, STBound, and MinIT on the linear relaxation of the three MIP formulations. • The impact of the above valid inequalities on the branch-and-bound algorithm used by Gurobi to solve the integer models. • How instance characteristics, determined by the instance generation parameters, affect the runtimes and optimality gaps of the MIP formulations. • How the CP approach compares with the best MIP model.

Impact of valid inequalities
In Section 5, we introduced five families of valid inequalities. All five of them are applicable to models M1 and M3, while only inequalities STBound and MinIT are valid for model M2.
We first measure the impact that each inequality, in isolation, has on the value of the continuous relaxation of the three models. To this end, we consider the value of the relaxation gap defined as Gap% = 100·(UB−LB)/UB where, for each instance, UB represents the best-known primal solution and LB is the objective value of the continuous relaxation. We present these results in Tables 2  to 4, where columns "T" report the time in seconds needed to solve the continuous relaxation and columns "#" list the number of valid inequalities added. Some columns "#" are omitted in Table 4 because the values for TaskInc, TaskSeq and CycElim are the same as in Table 2, while the values for STBound are the same as in Table 3. Values for MinIT are different between models M2 and M3 because the latter uses the dummy tray 0 while the former does not. Instances are grouped by two generation parameters, γ and δ, which have an important impact on the difficulty of the instances (see Section 7.1.3). Each row refers to the average over all instances sharing the  Table 2: Impact of valid inequalities on the strength of the continuous relaxation of model M1.
same parameter values for γ and δ. The last row reports, for each valid inequality, the percentage gap decrease and the root runtime variation, over all instances of the Synthetic set. The columns labelled with None, which report the results of the base formulations without valid inequalities, make clear that model M1 has the strongest relaxation, followed by M3. Model M2 exhibits large root gaps, which do not decrease adding valid inequalities.
The impact of valid inequalities on the other two models is positive. For M1, inequalities TaskInc have the strongest effect. Other inequalities only slightly decrease the root gaps. On the other hand, adding any type of inequality does not significantly increase the runtime at the root node. Note that we generate very few inequalities of family STBound; despite their number, they have a small but non-negligible impact on the root gap.
Valid inequalities have a larger effect on model M3, although not enough to make it competitive with M1. As for M1, inequalities TaskInc are effective; for model M3, however, CycElim are even more effective and achieve a root gap reduction of almost 0.6 percentage points. Inequalities MinIT slightly reduce the gaps but significantly increase the runtime, due to their high number.   Table 4: Impact of valid inequalities on the strength of the continuous relaxation of model M3.

Branch-and-bound performance
A tighter continuous relaxation does not necessarily correspond to a MIP which is easier to solve. Therefore, we use the results obtained in Section 7.1.1 as an indication, but still perform further experiments to determine the impact of valid inequalities on the overall branch-and-bound algorithm used to solve the three formulations. Because evaluating all possible subsets of inequalities to activate would be too time consuming, we add them incrementally. For each model, we start from the one which reduced the average percentage gap the most according to the results presented in Section 7.1.1, and proceed in decreasing order. For each model, we report results relative to two formulations. The first is the base model without valid inequalities, under header Base. The second is the formulation, among those we test following the above procedure, which gives the lowest average gap (see below for a precise definition of the gap). For model M2 the base formulation also gives the lowest gap. Table 5 reports the results of our experiments. Columns "Gap%" list the percentage optimality gap obtained with a time limit of one hour. The gap is defined as 100 · (UB − LB)/UB, where UB and LB are, respectively, the best primal and dual bounds obtained within the time limit. Columns "T" report the runtime in seconds, while columns "N" list the number of branch-and-bound nodes visited. Columns "Opt%" give the percentage of instances solved to optimality.
For model M1, we obtain the best results when adding inequalities TaskInc and MinIT; for model M3, when adding inequalities TaskInc and CycElim.
Model M1 clearly outperforms the other two formulations. When strengthened with the valid inequalities, in fact, M1 solves almost 99% of the Synthetic instances to optimality, with an average runtime of 72.021 seconds. Model M3 (with valid inequalities) also shows good performance, M1 Base Figure 6: Impact of instance generation parameters on the runtime of model M1 + TaskInc + MinIT.
solving more than 94% of the instances. Its average runtime and gap are worse than those of M1 over all subsets of instances and, therefore, we see no reason to prefer it to M1 even on a subset of the instances. Finally, formulation M2 is the weakest of the three. It starts with the weakest continuous relaxation at the root node and shows the highest gaps and the lowest number of closed instances at the end of the 1-hour time limit. Its average gaps, however, are always under 0.2%, showing that M2 suffers from a "tail effect": the dual bound increases quickly in the beginning, but then stops at a sub-optimal value. We also point out that additional experiments based on a dataset that was not generated with backbone instances (i.e., in which the instances were not necessarily feasible) showed that M2 was the fastest model to detect infeasibility. In practical applications, M2 could cover a complementary role to M1; for example, a decision-maker can run both models in parallel when it is not certain that a feasible elevator schedule exists.

Impact of instance generation parameters
We investigate which instance characteristics are associated with harder-to-solve problems. To this end, we use the best-performing model identified in Section 7.1.2, M1 + TaskInc + MinIT, and analyse how its runtime varies depending on parameters m, ω, δ, π gr , and γ. Figure 6 shows how average runtimes vary with the above parameters. Each plot corresponds to a different parameter and each box to a different value of that parameter. The average runtimes over all Synthetic instances are reported on the y axis, on a logarithmic scale. The horizontal lines inside the boxes and the numbers over them represent the median runtimes. Boxes span from the 1 st to the 3 rd quartiles; whiskers extend to the rest of the distributions, except for outliers marked with diamond flyers. A point is deemed an outlier if it lies more than 1.5 times the inter-quartile range outside the two central quartiles.
All instance generation parameters have a considerable impact on the solver runtime. Larger instances with more trays (higher values of m) are, predictably, harder to solve. Instances in Figure 7: Correlation matrix between instance generation parameters and key indicators related to the instances (the ratio of trays to shelves, m/n, and the average number of open time windows per instantĀ) and to model M1 + TaskInc (the number of rows NR, of columns NC, of non-zeroes NN, and the runtime in seconds). which jobs are "squeezed" in a small time horizon (small ω) and instances with longer shelf lives (smaller values of the greedy probability parameter π gr ) also correspond to longer runtimes. In general, instances with more overlap of task time windows are harder to solve, as confirmed by the large increase in runtimes when increasing the time-window half-width parameter γ. Increasing the number of tasks assigned to each tray (parameter δ) also makes the model more challenging.
To better study the relationship between instance and MIP model characteristics, we report in Figure 7 the correlation matrix between: • The instance generation parameters: m, ω, δ, π gr , and γ.
• Two indicators which measure how "dense" an instance is. The first is the number of trays per shelf (m/n). The second is the average number of open time windows per time instant, defined asĀ • Three indicators of the size of the resulting MIP: the number of columns (NC), the number of rows (NR, which includes both rows in the basic formulation of M1 and valid inequalities), and the number of non-zero entries in the constraint matrix (NN). • The runtime in seconds. Parameters m and δ have the largest impact on the size of the model, because they directly affect the size of set V , i.e., the number of y variables. IndicatorĀ also affects the size of the model because the more time windows overlap, the fewer tuples (t 1 , j 1 , t 2 , j 2 ) can be excluded from set V using the procedure described in Section 5.1. Larger models (higher NR, NC, NN) result in longer runtimes, while indicator m/n negatively correlates with runtimes. This is because, everything else being equal, in instances with larger m/n the same number of trays is packed in fewer shelves during the instance generation procedure. This is possible because tray shelf lives are short and non-overlapping, which are conditions which make fixing variables ∆ and y easier, and allow to generate more inequalities TaskInc.   Table 6 reports the results of this experiment. To compare the primal solutions from MIP and CP algorithms, columns "Gap%" refer to the percentage gap with the optimal solution of each instance (or the best-known solution for the only instance in the Synthetic set which we could not solve to optimality). We report the run-times in seconds in column "T". Column "Opt%" lists the percentage of instances for which the final solution returned by the algorithm was equal to the optimal solution. This does not mean that CP could prove that the solution was optimal: indeed, except for two small instances, CP could never prove optimality within the time limit. Even though CP is not competitive with model M1, its optimality gaps are small (around 1.5% on average). Further analysis on CP outputs showed that the best solution was usually found in the first five minutes. This indicates that, for very difficult instances, there could be an interest in running CP for a limited amount of time to find an initial solution and use it as a warm-start in the other models.

Results on the Realistic instances
This section reports results on the Realistic instances, which differ from the Synthetic ones due to their much larger time horizons and the higher number of tasks assigned to each tray. During preliminary experiments, we noticed that MIP-based models often struggled to find any primal solution. For this reason, during these experiments, we provided these models with an initial feasible solution obtained running the CP model M4 for 5 minutes (thus leaving 55 minutes as available runtime for the MIP models).
To account for the larger instance size, we increased the amount of allocated RAM from 4GB to 12GB. Even with this adjustment, though, models M2 and M3 ran out of memory on, respectively, 18 and 19 of the 50 instances. Therefore, we report results comparing models M1 + TaskInc + MinIT and M4, i.e., on the two models that solved all instances in the set. Moreover, model M2 proved that one of the instances we were provided with was infeasible. Thus, the following results refer to the remaining 49 instances. Table 7 reports the results of our experiments. Because the CP algorithm does not provide meaningful gaps and because we do not know the optimal solution to most of the Realistic instances, to obtain a uniform comparison we compute the values in columns "Gap%" as the percentage gap between the primal bound provided by each of the two algorithms and the dual bound provided by M1. Columns "T" report the average time in seconds. Each row aggregates results over all instances with the same number of shelves, whose number is reported in column "#". We note that the Realistic instances are hard to solve to optimality. Percentage gaps at the end of the time limit stay high, although it is not clear whether this is due to a low quality primal bound or a low quality dual bound (or both). We observe that for the largest instances (n ≥ 15), the primal bound found by the CP model was better than the one found by M1. Further analysis showed that for these instances, the CP model improved its best solution consistently until the end of the time limit, which is a major difference with respect to the Synthetic instances where the best incumbent was found within the first five minutes. We also note that the size of model M1 + was huge (around 35000 variables and 3.25 million constraints on average), mainly because of the valid inequalities. As we observed that most constraints were removed by the inner preprocessing of Gurobi, we tried the simplest version of M1 (i.e., without any valid inequalities), but the results were significantly worse.

Conclusions and future research
We introduced the Vertical Farming Elevator Energy Minimisation Problem (VFEEMP), a realworld problem arising from the operations of autonomous vertical farms. We showed that its decision version is NP-complete and presented three MIP formulations together with a set of valid inequalities, and a CP model for the problem. We also introduced two sets of benchmark instances, one derived from real-life data, and the other to determine the instance characteristics that make the problem harder to solve. We empirically evaluated the performance of the models and the impact of the valid inequalities and showed that (i) MIP model M1 with valid inequalities TaskInc and MinIT obtained the best performance on synthetic and small-size realistic instances, and (ii) CP model M4 helps finding good quality initial solutions for small-size synthetic instances and was superior to M1 for large-size realistic instances. We also identified parameters that make instances harder to solve for our models: a large number of trays, a large number of tasks per tray, and large time windows. Indeed, the two first characteristics directly influence the number of variables and constraints in the models, while the latter increase the number of possible quadruple (t 1 , j 1 , t 2 , j 2 ) that are valid indices for variables y.
Future work shall focus on the development of specialised heuristics for the problem. It would also be interesting to study the development of decomposition approaches. For example, disregarding time window violations in the master problem and checking the validity of the resulting schedule in the slave problem. Another possibly fruitful research avenue is time decomposition. For example, in Figure 5, solving the instance first for trays 1,2,3,5, and the first tasks of tray 4, and then tackling the remaining tasks of tray 4 and tray 6. However, since the windows are relatives, such decomposition would not necessarily reach an optimal solution.
Furthermore, other relevant real-world variants of the problem could be studied: in a generalisation of our problem, the planner could be allowed to decide the shelf associated with each tray. In other words, the control system of the VF tower would have to jointly assign trays to shelves and plan the movements of the elevator. In a further generalisation, we could allow the elevator to move a tray from one shelf to another (empty) shelf at any point during the growth period of the tray. Even though this operation has some energy costs, it can potentially bring more savings on the long run, for example, when bringing an isolated tray closer to the rest of the occupied shelves.