Solving the Traveling Telescope Problem with Mixed-integer Linear Programming

The size and complexity of modern astronomical surveys has grown to the point where, in many cases, traditional human scheduling of observations are tedious at best and impractical at worst. Automated scheduling algorithms present an opportunity to save human effort and increase scientific productivity. A common scheduling challenge involves determining the optimal ordering of a set of targets over a night subject to timing constraints and time-dependent slew overheads. We present a solution to the “traveling telescope problem” that uses mixed-integer linear programming. This algorithm is fast enough to enable dynamic schedule generation in many astronomical contexts. It can determine the optimal solution for 100 observations within 10 minutes on a modern workstation, reducing slew overheads by a factor of 5 compared to random ordering. We also provide a heuristic method that can return a near-optimal solution at significantly reduced computational cost. As a case study, we explore our algorithm’s suitability to automatic schedule generation for Doppler planet searches.


INTRODUCTION
Maximizing the scientific yield of expensive and often oversubscribed astronomical instrumentation requires meticulous planning of each night's observations.However, determining the optimal (or even near-optimal) sequence of observations is challenging and time consuming task.Schedulers must incorporate temporal accessibility windows while factoring in slew and acquisition overheads that can themselves be time-variable.In this paper, we refer to the task of determining the optimal ordering of a set of observations by a telescope as the 'Traveling Telescope Problem' or 'TTP' given its similarities to the 'Traveling Salesman Problem' or 'TSP'.
The scientific benefits of intelligently sequenced observations can be significant, especially for programs with many targets spread over the entire sky.As an example, the Doppler planet searches at the Keck-I telescope observe up to 100 targets per night (Howard et al. 2010).As we show below, a quasi-random sequence of 100 targets requires over 3 hours of slew during a 10 hour night while an optimized sequence can reduce this to 0.5 hours.
Automated solutions to the TTP offer a number of opportunities.As is the case for the TSP, a skilled human scheduler can generate a observation sequence that significantly outperforms a random ordering.How-ever, such script generation takes significant human effort that could be devoted elsewhere.In addition, a number of ongoing and forthcoming surveys are or will be scheduled automatically.The Zwicky Transient Facility (ZTF; Bellm et al. 2019a) and the Legacy Survey of Space and Time (LSST; Ivezić et al. 2019) are two examples.Automated solutions to the TTP are necessary for automated surveys.
While a rich literature exists on the TSP, standard solutions do not directly transfer to the TTP for several reasons.First, targets are often only accessible for a fraction of the night either due to their position on the sky or scientific need for time-critical observations.Thus, a large fraction of the N !target sequences are infeasible.Second, slew time between targets is itself a function of time.As an example, Figure 1 shows the trajectory of two targets above the Keck-I telescope atop Maunakea.Two targets cross the meridian north and south of zenith, respectively.Like most large telescopes, Keck-I moves in the altitude and azimuth directions and target slews are dominated by differences in azimuth.Figure 1 shows the variation in slew time over the course of the night, which grows as the targets first approach the meridian and then straddle the telescopes cable wrap limits.
Previous efforts have addressed the TTP under a number of simplifying assumptions.The ZTF scheduler divides the night into short intervals and solves a standard TSP while treating the set of accessible targets and target-to-target slew times as constant with in the interval (Bellm et al. 2019b).The James Webb Space Telescope's scheduling architecture (Giuliano & Johnston 2008), which is built upon the Hubble Space Telescope's SPIKE software (Johnston & Miller 1994), also treats target-to-target slew times as constant.At present, we are not aware of any global solutions to the TTP that capture both variable accessibility windows and slew overheads.
This paper is organized as follows.We describe the TTP in Section 2 and present a formulation using Mixed-Integer Linear Programming that can be solved to global optimally range of problem sizes.In Section 3 we conduct a suite numerical experiments to determine the performance and computational cost of our algorithm over various problem sizes.Section 4 discusses the limitations of our global approach and the potential of heuristic solutions.We conclude in Section 5.

Problem Setup
For a given set of targets and an observing interval, we seek the tour that completes all exposures in the shortest possible time.Targets must be accessible for at least part of the observing duration.The slew time between any pair of targets must be computed in advance, but the slew time may itself be a function of time.The formulation presented below was inspired by Sun et al. (2018) who developed a framework to optimize the profitability of parcel pickup and delivery with variable time windows and travel times.

Slot Framework
Following the TSP literature, we refer to targets to be traversed in the TTP as nodes since that work emphasizes that the solution is a directed graph connecting all targets.With a list of N targets to be scheduled, we define the total set of nodes to be {0, 1, ..., N, N + 1}.The nodes 0 and N +1 are the anchoring start and end nodes; their location is arbitrary, i.e. not associated with a celestial source.Their purpose is explained in Section 2.4.Each node i has an associated exposure time τ exp i , and accessibility window [t e i , t ℓ i ].The values t e i and t ℓ i are the earliest and latest times the tour can depart node i, i.e. the time the observation concludes (see Figure 2).We summarize the symbols from the main body of this text in Appendix A.
Next, we break the scheduling interval into M subintervals or 'slots', within which travel time is treated as constant.M is a free parameter and impacts the computational load (see Section 3).The slots may have uneven durations if desired.The boundaries of the slots are [w m , w m+1 ] where m indexes each slot.
We then construct the slew tensor τ slew ijm that specifies the travel time between any two nodes during every slot (see Figure 1).This depends on the telescope slew speed in altitude and azimuth directions as well as cable-wrap considerations.For a concrete example, τ slew 4,6,2 specifies the computed travel time between the i = 4 node and j = 6 node during the bounds of the slot m = 2.

Decision Variables
Our MILP formulation involves both binary and continuous variables.The following variables trace the flow through the nodes and the times of node departures: • X ijm is a binary variable equal to 1 if the tour traverses the arc from node i to node j during slot m and 0 otherwise.
• Y i is a binary variable equal to 1 if the node i is entered at some point during the tour and 0 otherwise.
• t i is a continuous variable denoting the departure time from node i.
• t ijm is a continuous variable and equal to t i if the tour departs the node i toward the node j during the time slot m, 0 otherwise.
The non-zero elements of X ijm describe the tour.The tensor dot product of X ijm and τ slew is equal to the total travel time.

Constraints
Next, we introduce the following constraints: Constraint 1. Tour must depart the starting node.We require that flow be non-zero from 0 to some arbitrary target node j during some slot m in time to begin the tour.
Constraint 2: Tour must conclude at the ending node.Similarly, we anchor the end of the tour with the end node N + 1.We must traverse the arc (i, N + 1) from some arbitrary target node i.
Constraint 3: Y must indicate the visitation of a node.In the simplest possible case, we have a trivial solution traversing from 0 to N + 1, stopping at a single target node along the way.The objective (see Section 2.5) will reward the visitation of additional nodes between these two constrained anchors.This requires first defining the variable Y to track whether nodes are being visited.Y i will be activated if the tour flows from the node i during any slot in time, into any subsequent node j.
Constraint 4: A slew into any target node must be accompanied by a slew away from that node.
Excluding the starting and ending nodes, we require that the slew into any target node k must be accompanied by a slew out of k.
(5) If X ijm is 1 for any arc (i, k) in the sequence, it will also hold 1 for some (k, j) arc at an arbitrary time.If the left term is 0 (the tour never enters node k) then the tour will never traverse a departing arc originating at k.The chronology of these events will next be enforced using the time variable t.Constraint 5: Define the selection variable t ijm .Now we enforce the proper constraints to define the selection time variable t ijm relative to our generic time variable t i .By summing over all potential destinations j and time slots m, we recover the value stored in t i .
The second time variable t ijm is critical for the following two constraints.It behaves like the product of t i and X ijm , but works within a linear program.Constraint 6: Departure time from a node is greater than the departure from the previous node plus the slew and exposure time.Say we begin to traverse from one node to another on the arc (i, j) We have exaggerated slew times to clearly illustrate the differences between the three tours.
In the top sub-optimal tour, we include a horizontal timeline for targets A-D.The windows of accessibility are shown as vertical ticks.For example, the earliest the telescope can complete observations of target A is labeled with t e A while t l A , denotes the latest time.Note that t e A depends on both the rise time and exposure duration, and t ℓ i is the either the set time or the end of the observing interval, whichever comes first.This scheduler uses a greedy approach; i.e., once the observation A of completes, the telescope immediately slews to the unobserved target with the shortest slew time.The tour is feasible, but the entire observing duration is required in this example.The middle optimal tour is scheduled with a global optimization.In our example, the slew between A and D is so long that it is advantageous for the telescope to wait until target B rises to observe it, before proceeding to D. This is the optimal tour to observe targets A-D.Finally the bottom plot shows the inclusion of a fifth target E to illustrate an infeasible tour.There is not enough time to observe all targets.
within the bounds of the slot m.Before the telescope may depart from the next node j, a minimum amount of time must pass, equal to the slew experienced on the journey from i to j plus the exposure time at the new node, i.e. τ slew ijm + τ exp j .The minimum time value of our departure from j is the time that the previous node i was departed t ijm plus this minimum 'passing time'.The inclusion of the variable X ijm will ensure only the proper values of i and m are considered for the journey to the new node j.
Constraint 7: Ensure departure times are consistent with slot bounds.We must enforce constraints on the departure times t ijm using bounds of the time windows defined in w.If we exit the node i during the time slot m, then the departure time must be within the bounds of the slot window.
Constraint 8: Departure times must respect node accessibility.Finally, we enforce the time window constraints on the departure times t i for all the visited target nodes. 2.5.Objective We seek to maximize the the number of scheduled exposures while minimizing total slew time Here, C is a small constant such that the slew penalisation (second term) is always less than unity.Notice that our anchor nodes do not contribute to the total slews with this summation convention, and are therefore used only for constraining the flow of the tour (a physical location need not be set, and the values in the first row and column of τ slew ijm are set to 0).The objective function does not require all nodes be visited.In similar TSP literature such as Sun et al. (2018), each target node may be assigned a scalar priority p j included as a coefficient in the left summation term in Equation 10.In such a formulation, lowest priority targets are preferentially excluded when total completion is infeasible.We note that global optimality becomes less intuitive when targets have different numerical priorities.
In the TTP, the optimization will work to remove the targets that contribute most to the total slews in every case.The objective is clear: observe as many targets as possible as quickly as possible.We note that the formulation above describes a simple TTP where targets are observed once and no additional constraints are placed on the timing between observations.Appendix B describes small modifications to the algorithm presented above that can accommodate such constraints.

Numerical Experiments
We evaluated our algorithm's ability to solve the TTP through a suite of simulated target lists.We explored problem complexity along the following three axes: number of targets N , number of time slots M , and duration of the scheduling interval D. We designed TTPs in the following manner: 1. We specified a random calendar night at Keck Observatory.
2. We selected an observing duration D.
3. We selected N stars from the California Legacy Survey (CLS; Rosenthal et al. 2021), a collection of 719 nearby stars that have been observed from Keck observatory for several decades as part of an extrasolar planet search.These stars comprise a good TTP test set because they are nearly uniformly distributed on the sky with declination δ ≳ −40 deg and are thus observable for ≳ 7 months per year.When sampling the targets, we required they be accessible for more than 50% of duration D. We removed a few close binaries that have negliable slew speeds.
4. We modeled the slew rate of the Keck telescope as 1 deg per second in both altitude and azimuth directions.
5. We assigned uniform exposure times (in minutes) to all targets according to: Setting the exposure times this way would completely fill the observing duration given a random ordering of targets with an average slew time of 120 sec or 2 min.For our experiments, slews in azimuth (which ranges from 0 to 360 deg) dominate over those in altitude (which ranges from 0 to 90 deg).The average distance between two randomly selected values between 0 and 360 is 120.We expect all solutions to the TTP to be significantly more efficient.Thus our experiments are feasible by construction and we report the improvement relative to this random ordering.
6.We selected M uniform slots within D.
7. We calculated the distance tensor τ slew ijm at the midpoint of each slot.
With the experiment specified, we solved the MILP described in Section 2 with one additional constraint.Constraint 9: All target nodes must be visited.
The problem generation process described above in general results in TTP instances in which it is possible to observe all N targets.Given this, we modify the objective function described in equation ( 10) to focus on minimizing slews only, resulting in the following new objective function.

Min
In our results in Section 3.2, we use TTP-Global to refer to the formulation defined in Section 2 with the objective function in equation ( 13) and constraint (12).
Before describing our grid-based exploration of problem complexity, we show one example solution to the TTP in Figure 3.We compared it with a simple heuristic solution for a 50 target observing sequence conducted over a half night to emulate a human generated script.In this heuristic, targets are observed in the order of their set times, i.e. the earliest setting target is observed first.Figure 3 shows graphically the slew overheads that TTP-Global eliminates through a more efficient ordering.

Computational Results
We test our formulation using different combinations of the number of targets N , the number of slots M and the duration/scheduling interval D. For the scheduling interval D, we consider quarter nights, half nights and full nights.For quarter nights, we vary N in {5, 10, 25}; for half nights, we vary N in {5, 10, 25, 50}; and for full nights, we vary N in {5, 10, 25, 50, 100}.For each combination of D and N , we vary M in {1, 3, 10}.
For each combination of N and D, we consider 10 randomly generated sets of targets, and consider the three different values of M , giving rise to a total of (3 + 4 + 5) × 3 × 10 = 360 problem instances.We solved TTP-Global using Gurobi version 10.0.1, a state-of-theart optimization suite that solves MILP problems using the branch-and-bound algorithm (Gurobi Optimization, LLC 2023).We used the Python programming language to generate the input data for TTP-Global and to formulate TTP-Global using the Gurobi Python API.We conducted our suite of numerical experiments on Amazon Elastic Compute Cloud (EC2), on a single instance of type m6a.48xlarge (AMD EPYC 7R13 processor, with 192 virtual CPUs and 768 GB of memory).For each experiment, we allocated 8 virtual CPUs and limited computation time to 1800 seconds.
For a few points of reference, a TTP with N = 25 targets and M = 1 involved an MILP with 1643 rows (constraints) and 1513 columns (variables), and the M = 10 case had 14765 rows and 14635 columns.A TTP with N = 100 targets had 21518 rows, 21013 columns for M = 1, and 208790 rows, 208285 columns for M = 10.
For each experiment, we recorded the following information: • SlewTime τ : total slew time of the final schedule obtained from Gurobi, calculated using the discretized tensor τ slew ijm .• RelRed τ : relative reduction in slew time of the final schedule compared to the randomly ordered value of 2N ; mathematically, it is defined as: • SlewTime real : real slew time of the final schedule, calculated based on the t i departure time values of the schedule.
• RelRed real : the analog of RelRed real for the real slew time.
• Runtime: computation time • Whether a provably optimal solution was found.
• Whether a feasible solution was found.
Table 1 summarizes these statistics for the ten experiments conducted at each combination of N , M and D. We report the number of experiments where Gurobi found a feasible solution, NumFeas, and a provably optimal solution NumOptimal.For the feasible set, we report the average values of SlewTime τ , RelRed τ , SlewTime real , and RelRed real .for each combination of N , M and D, where the average is taken over the NumFeas instances for which a feasible solution was found.For example, for (Half, 50, 3), NumFeas is 8, indicating that Gurobi found a feasible schedule in only 8 out of the 10 instances; consequently, the value of 12.7 for SlewTime τ is the average slew time over those 8 feasible schedules.We show the average Runtime and RelRed real values for different problem sizes in Figure 4.

DISCUSSION
We may draw a number of conclusions about the suitability of TTP-Global for the TTP from Table 1 and Figure 4. Gurobi generally found a feasible schedule when the number of targets N ≤ 25.For N ≥ 50 targets, our ability to find a feasible solution depended sensitively on the number of slots.Gurobi found a feasible schedule for all ten instances when M = 1, for some when M = 3, and for none when M = 10.
Not surprisingly, runtime was a strong function of N and M as can be clearly seen in Figure 4.For example, Gurobi found optimal solution for all (D, N, M ) = (Full, 25, 1) experiments with an average Runtime of 4.7 seconds.In contrast, the (Full, 25, 10) experiments found no optimal solutions in the 1800 second time limit.For the largest (Full, 100, 3) experiments, not even a feasible solution was found in the time limit.We find that the M = 1 case scales well into the realm of N = 100, solving even faster than the N = 50 case.While this goes against our initial intuition, we suspect this behavior is the result of parameter tuning by Gurobi to accommodate larger models.
When we do find a feasible schedule, it is significantly more efficient relative to the 2N baseline.For example, for the (D, N, M ) = (Half, 50, 3) experiments, the average SlewTime real is 22.4 minutes compared to the 2N = 100 minute benchmark, a reduction of RelRed real = 77.6%.For the (Full, 100, 1) set of instances, average SlewTime real is 36.1 minutes, which is a reduction of 82.0% relative to the 2N = 200 minute benchmark.We note that in all cases, SlewTime τ is less than SlewTime real .For example, for the same D = Full, N = 100, M = 1 set of instances, SlewTime τ is a mere 16.4 minutes, which is a reduction of RelRed τ = 91.8%relative to the 200 minute worst-case value.The difference between RelRed real and RelRed τ stems from our piece-wise slew approximation, i.e. the choice of M for each model.
In the quarter night models, we find little improvement from using higher values of M , since the slews vary minimally with respect to time.For the longer durations, we see some benefit for higher M in the simpler N < 25 cases.For N ≥ 25, the higher M models be-  1. Computational results of TTP-Global for different values of D, N and M .Note: "-" indicates that the metric could not be calculated due to Gurobi not being able to obtain a feasible schedule for any of the ten instances.For D = Full, N = 5 and D = Full, N = 10, one instance was determined to be infeasible.The goal of our experiment is to devise target lists and exposure times that are feasible in the limit of large N, however for small N this is not strictly guaranteed.
come too complex to be solved to optimality in the time limit (leaving better slews in question), but the M = 1 models continue to demonstrate extremely efficient slews despite the higher expected variability.
In most cases, increasing the value of M did not demonstrate a significant advantage over the static models.For the high N models, the optimizer was often not able to construct a feasible solution for M > 1.For scheduling full nights of observations, the M = 1 case demonstrates dramatic improvement in slew times by up to a factor of 5, and increasing M provides more computational complexity than can be handled by our global algorithm in the time limit.For most cases, we slews by around a factor of 5.For small N , there is some benefit to using M > 1.For larger models, the risk of long slews in the M = 1 case has little impact on RelRedreal.
would recommend M = 1 be the standard due to its exceptionally fast runtime.
In our numerical experiments we attempted to solve TTP-Global to a provably optimal solution by branchand-bound.For large numbers of targets or finely discretized slew tensors this approach may not be computationally tractable.We note that many heuristic solutions to the TSP have been developed that achieve near-optimal solutions.We suspect that analogous highquality heuristic solutions exist for the TTP, which may be equipped to solve the M > 1 models for much higher N .We develop one such procedure in the Appendix for comparison with our global solution.

CONCLUSIONS
In this work, we addressed the challenge of determining the most efficient ordering of a set of exposures at a telescope.We formulated the Traveling Telescope Problem (TTP) as a mixed-integer linear program TTP-Global and solved it using a standard commercial optimizer for problem sizes as large as 100 targets in ∼ 10 minutes using modest computational resources.The speed of TTP-Global means it can be run dynamically throughout the night and respond to real-time changes in target accessibility from weather.Further work is needed to develop algorithms that can solve the TTP for substantially larger sets of targets or substantially finer time-resolution in slew overheads.Local searches initialized with an heuristic solution may prove fruitful.We hope that algorithms like the one described here can assist or automate scheduling, save human effort and increase the scientific productivity of astronomical surveys.
L.H. & E.A.P. acknowledge support from the following sources: Heising-Simons Foundation Grant #2022-3832.V.V.M. acknowledges support from the UCLA Anderson School of Management.We are grateful for enlightening conversations with Eric Bellm.
The first constraint ensures that if both observations take place, the directed tour must pass directly from i to i ′ or vice versa.The second constraint ensures both observations or neither observation occur.

B.2. Intra-Night Spacing Requirements
One may wish to enforce a minimum interval between two exposures.A common example occurs in time-series monitoring where multiple observations of the same target occur during the same night, subject to a minimum separation.We accomplish this by letting N correspond to the total number of exposures to be collected across all targets.For simplicity let us assign exposure indices such that exposures of the same target occur consecutively in the total exposure list {1, . . ., N }.That is, if the target requires n exp individual exposures, the node indices {κ, κ + 1, . . ., κ + n exp − 1} correspond to that target for some value κ.
With the repeat observations specified, we enforce a minimum interval via the following constraint: Subsequent exposures of a given target must not end until at least τ sep has passed since the previous exposure ended.
For example, if the linked exposures of a target have index i = 5 and 6, exposure 6 may not end until at least τ sep has passed since exposure 5 ended.

C. LOCAL SEARCH HEURISTIC
Here, we develop an alternative formulation to the TTP which uses a local search heuristic.We refer to this as TTP-LS to distinguish it from TTP-Global which searches the global solution space.

C.1. Algorithm Description
In the TTP, there are two sets of decisions that need to be made simultaneously.One is the sequence in which the targets will be visited.For example, with N = 5, we have to choose between 5 → 1 → 4 → 2 → 3, 1 → 3 → 2 → 4 → 5, and so on.The other set of decisions involves the timing slews.This involves deciding a (continuous) time t i within the observing duration D to slew and the corresponding slot m.Even with a fixed sequence exposures, this is a non-trival task.As a result, making both sets of decisions under the umbrella of a single MILP formulation is computationally demanding.
This suggests an alternate approach to the TTP that decouples the sequence decision from the timing decision.Suppose that the sequence of targets is fixed.When should the telescope slew in order to minimize slew times?Let σ denote the sequence of targets, which is a bijective function σ : {1, . . ., N } → {1, . . ., N }, and let the minimum total slew time be denoted by the function F , so that F (σ) is the minimum total slew time that one would obtain from following the sequence σ.Let Σ denote the set of all such sequences.For example, for N = 5 targets and the sequence 5 The TTP can then be abstractly formulated as min σ∈Σ which is an optimization problem over sequences in Σ.As written, this is not a problem that can be readily provided to any commercial solver, but because of the discrete nature of Σ, it can potentially be solved using local search.Let z ∈ {1, . . ., N } and z ′ ∈ {1, . . ., z − 1, z + 1, . . ., N } be positions in the sequence, and let σ z↔z ′ denote the sequence obtained by swapping the targets in positions z and z ′ ; that is, σ z↔z ′ is the unique sequence such that Let N z (σ) denote the set of neighboring sequences of σ obtained by swapping the target at position z with a target at any other position: With this definition, our local search algorithm can be formally described as Algorithm 1.
In words, we begin from some initial sequence σ.We use U to denote the set of sequence positions which have not yet tried to modify.As long as there is at least one sequence position we have not tried to change, we pick a sequence position z, and calculate the best neighboring sequence σ * obtained by swapping the target at position z with the target at any other position.If the objective value F ′ of the best neighboring sequence improves on the current objective value F (σ), we replace σ with σ * , and we reset U to be the set of all positions.Otherwise, if we do not make improvement in an iteration of the while loop, then U will be reduced by one member.If N such iterations occur, then U will be empty, and we will have ascertained that there is no neighboring sequence we can move to in order to reduce the objective value; in other words, σ is a locally optimal sequence.We note that this heuristic is similar to the 2-OPT heuristic (Croes 1958) for the classical TSP problem, which involves eliminating two edges in a end if 10: end while 11: return Locally optimal sequence σ Algorithm 1: Pseudocode of local search procedure.
TSP tour and reconnecting the tour so that the edges are swapped.The main difference comes from the function F , which calculates the minimum total slew time when one assigns the targets in the sequence to slots optimally.
Before this algorithm can be deployed, we must specify how to compute F (σ).The function value F (σ) for a fixed sequence σ can be calculated by solving a MILP.While this MILP shares some similarities with the TTP problem described in Section 2, it is a simpler because target sequence is fixed and "baked in" to the optimization problem.This smaller MILP is a subroutine in Algorithm 1.We provide further details on this integer program in Section C.2.
We must also consider sequences of targets that are infeasible.In some cases, for a fixed sequence σ of targets, it may not be possible to make the timing decisions and the slot decisions in way that respects the accessibility windows and slot bounds.In such cases, the MILP that defines the function F (•) will be infeasible.We can extend the definition of F (•) so that F (σ) = +∞ if the corresponding MILP is infeasible for sequence σ; since Algorithm 1 is always choosing the neighboring sequence with the lowest value of F , this ensures that Algorithm 1 will never replace the current sequence with one that is infeasible.
However, even with this fix, one problem that still remains is if the initial sequence σ and all neighboring sequences of that initial sequence are infeasible.In this case, Algorithm 1 will not return a feasible sequence, as it will simply terminate with the current sequence.This is a serious issue, because it is not straightforward to identify a sequence of targets for which the TTP problem constraints can be perfectly satisfied.In Section C.3, we present a feasibility heuristic for identifying such a sequence.

C.2. Calculating Minimum Total Slew Time for a Fixed Sequence of Targets
As noted in the previous section, a key component of Algorithm 1 is the function F , which maps a sequence σ to a minimum slew time F (σ).We use z to denote the index of a position in this sequence σ.For targets, z will range in {1, 2, . . ., N }.With a slight abuse of notation, we will use z = 0 to denote the start node of the telescope, and assume that σ(0) = 0; similarly σ(N + 1) = N + 1 at the final node.Thus, z can take any value in {0, 1, . . ., N + 1}.
Let Y at z,m be a binary decision variable that is 1 if the telescope departs the target at position z ∈ {0, 1, . . ., N + 1} in the sequence in slot m, and 0 otherwise.Let Y by z,m be a binary decision variable that is 1 if the telescope reaches slot m by position z ∈ {0, 1, . . ., N + 1} in the sequence and 0 otherwise.Let t z denote the departure time of the telescope from the node at position z.The function F is obtained by solving the following MILP: subject to Y by 0,0 = 1, (C4b) t z ≥ t e σ(z) , ∀ z = 0, 1, . . ., N + 1, (C4j) Y by z,m ∈ {0, 1}, ∀ z = 0, 1, . . ., N + 1, m = 0, 1, . . ., M − 1, (C4l) In order of appearance, the constraints have the following meaning.Constraint (C4c) requires that if we have reached slot m by position z, then it must be the case that we have reached slot m by position z + 1. Constraint (C4d) requires that if we have reached slot m+1 by position z, then we must have reached slot m by position z.Constraints (C4e) and (C4f) link the Y by and Y at variables; constraint (C4e) means that we are in slot m at position z if and only if we have reached slot m by position z (Y by z,m = 1) and have not reached slot m + 1 by position z (Y by z,m+1 = 0).Constraint (C4f) similarly requires that we are in slot M − 1 at position z if and only if we have reached slot M − 1 by position z.Constraint (C4g) requires that the departure time from the target in position z is at least the slew time that is realized departing from the target in slot z − 1 plus the exposure time of the target in position z.Constraints (C4h) and (C4i) ensure that the departure time of each position z is within the lower and upper bounds of that position's assigned slot, while constraints (C4j) and (C4k) ensure that the departure time from each position z ∈ {1, . . ., N } is within the rise and set times for the target in that position (t e σ(z) and t ℓ σ(z) respectively).The last two constraints enforce that the Y at and Y by variables are binary.
The MILP problem (C4) is essentially the TTP problem of Section 2, restricted to a particular sequence σ.Essentially, this formulation decides when each target's departure time will be and to what slots the different positions in the sequence will be allocated.Importantly, the sequence of targets is not a decision variable as it is in the original TTP model; it is a fixed input that is provided by the user.
Many of the constraints are direct analogs of constraints that appear in the TTP-Global formulation.For example, (C4j) and (C4k) model the rise and set time constraints for each target in each position, mirroring constraint (9) of the TTP MILP.As another example, (C4i) and (C4h) model the lower and upper bounds of each slot, similarly to constraint (8).Note that because the sequence of targets is fixed, many of the constraints from the full TTP MILP can be simplified, and other decisions, such as the departure times and in which slot each target is being departed from, can be expressed more efficiently using different decision variables.Specifically, the decision of which slot each position is assigned to is captured by the Y by z,m decision variables.Here, we remark that these variables are an example of the incremental encoding technique in integer programming, which enhances the efficiency of branching in the branch-and-bound algorithm that is the cornerstone of integer programming solvers.We refer interested readers to the review paper of Vielma (2015), and to Bertsimas et al. (2011), Bertsimas et al. (2019) and Mišić (2020) for examples of applications of this technique in air traffic control, vehicle routing and optimization over trained machine learning models.
As a result, problem (C4) is much easier to solve than the full TTP MILP.We found that that Gurobi could determine the exact optimal solution to this problem in under a second with a single thread.

C.3. Feasibility Algorithm
The local search approach described above may fail if it initialized at an infeasible sequence whose neighbors are also all infeasible.In order for Algorithm 1 to return a sequence that can be implemented, the initial sequence must be one for which the minimum slew problem (C4) is feasible.Given the combinatorical complexity of the TTP, it is unlikely that one would randomly select a target sequence that would result in problem (C4) being feasible.
Thus motivated, we present in this section an algorithm that, starting from any sequence, seeks to return a feasible sequence.We note that this algorithm is a heuristic, and is not guaranteed to succeed.Nevertheless, our numerical results in Section C.4 indicate that this heuristic is generally very effective.
At a very high level, the algorithm we will propose resembles our local search procedure, Algorithm 1, in that it starts from a sequence σ and makes moves to neighboring sequences.The key difference is the objective function that is used.Instead of using the function F , the feasibility algorithm first seeks to locally optimize a function G 1 , followed by a function G 2 .The function G 1 (σ) measures, for the sequence σ, the smallest violation of the slot window constraints (C4h) and (C4i) that can be attained when we choose the departure times and the slots to which each target is assigned to.This violation is a nonnegative quantity; a positive value implies that we are unable to satisfy all of the constraints, i.e., at least one constraint in constraint sets (C4h) and (C4i) is violated.A value of zero implies that all of the constraints in the two constraint sets are satisfied.The function G 2 (σ) similarly measures the smallest possible violation of the visibility constraints (C4j) and (C4k) when we choose the departure times and the slots.Again, a positive value implies that at least one constraint in the constraint sets (C4j) and ( C4k) is violated, while a value of zero implies we can satisfy all of the constraints defined by ( C4j) and (C4k).
G 1 is defined by the following MILP: minimize t z ≥ t e σ(z) − ϵ e z , ∀ z = 0, 1, . . ., N + 1, (C5k) Observe that this integer program is similar to the minimum slew problem (C4), except that we now allow for violations of the constraints (C4h), (C4i), ( C4j) and (C4k).The main modifications are as follows.First, observe that in addition to the decision variables of problem (C4), problem (C5) includes the decision variables ϵ e z , ϵ ℓ z , ϵ slot,ℓ z , ϵ slot,e z , which measure how much the rise, set, slot upper bound and slot lower bound constraints can be violated in time.
Second, observe that constraints (C5i) -(C5l) resemble constraints (C4h) -(C4k), except that the new ϵ e z , ϵ ℓ z , ϵ slot,ℓ z , ϵ slot,e z appear.These new decision variables are bounded from below by zero and unbounded from above, so they effectively allow the optimizer to choose to not satisfy these constraints.For example, in the constraint t z ≤ t ℓ σ(z) + ϵ ℓ z , for whatever value of t z we choose, we can always make the constraint satisfied by setting ϵ ℓ z to be equal to any value greater than max{0, t z − t ℓ σ(z) }; as a concrete example, if t z = 120 and t ℓ σ(z) = 80, then any ϵ ℓ z ≥ max{120 − 80, 0} = 40 will satisfy the constraint.Lastly, observe that the objective function is equal to the sum of the ϵ slot,ℓ z and ϵ slot,e z variables.Thus, the optimizer seeks to find the assignments of sequence positions to slots and the departure times so as to minimize how much the slot lower and upper bound constraints from problem (C4) are violated.Observe that if the optimal objective value of problem (C5) is zero, then we have found a partially feasible solution to problem (C4) that satisfies all of the constraints, and in particular constraints (C4h) and (C4i), with the possible exception of the rise and set time constraints (C4j) and (C4k).Importantly, note that no matter what sequence σ one chooses, problem (C5) is always feasible.
We now define the function G 2 .The function G 2 is defined as the objective value of the following integer program, which is With problems (C5) and (C6) defined, we can now define the feasibility algorithm, which we provide in Algorithm 2. This algorithm works by first performing local search using the function G 1 .If the local optimum is such that the value of G 1 is positive, then the algorithm terminates and returns that the problem is infeasible.Otherwise, if the value of G 1 is zero, then we continue to the next phase, in which we perform local search using the function G 2 .If the local optimum of G 2 is positive, then the algorithm again terminates and returns that the problem is infeasible.Otherwise, if the value of G 2 is zero, then we have identified a feasible sequence.Note that Algorithm 2 is a heuristic and does not provably verify that problem (C4) is infeasible.If it returns "Problem is infeasible", it may not be the case that the minimum slew problem is actually infeasible.

C.4. Computational results for local search heuristic
We now present our results on our heuristic approach described above.We tested our approach on the same collection of 360 experiments described in Section 3.2, and compute the same result metrics.We tested two variants of our local search procedure: • TTP-LS-1: Here, we execute our overall algorithm from a single random starting point, which we obtain by drawing a sequence σ uniformly at random from all possible N !sequences.
• TTP-LS-10: In the second variant, we execute our overall algorithm from ten randomly generating starting points, each of which is a uniformly randomly generated sequence, and retain the best solution obtained over the ten repetitions.
With both TTP-LS-1 and TTP-LS-10, we impose a time limit of 600 seconds on the total run time.In the most extreme case, TTP-LS-1 will require 600 seconds, while TTP-LS-10 will require 600 × 10 = 6000 seconds.In both variants, the functions G 1 and G 2 (see Appendix C.3) and F (see Appendix C.2) are computed by solving the corresponding MILPs using Gurobi with a single thread.We again implement our procedure in Python and run our experiments on the same Amazon EC2 instance described in Section 3.2.Table 3 summarizes the results for TTP-LS-1 and Figure 5 shows run time and slew efficiency for different problem sizes.TTP-LS-1 exhibits favorable performance in terms of computation time; in most cases, TTP-LS-1 terminates with a locally optimal solution within 600 seconds.The only exception is the (D, N, M ) = (Full, 100, 10) set of instances.Note that TTP-LS-1 resulted in a feasible schedule in all but seven of the 360 instances; importantly, TTP-LS-1 finds a feasible schedule in all of the instances for parameter combinations for which the TTP-Global MILP fails (for example, for (Full, 100, 10), TTP-LS-1 produces a feasible schedule in all ten instances in 600 seconds, whereas TTP-Global fails to find a feasible schedule in all ten instances with 1800 seconds of computation.Of those seven instances in which TTP-LS-1 did not find a feasible schedule, six are the same instances which were determined to be infeasible by TTP-Global, and in one instance, the feasibility procedure (Algorithm 2) failed to identify a feasible solution, despite the fact that the instance does admit a feasible solution based on running TTP-Global.Lastly, in terms of solution quality, the total slew time, as measured by SlewTime τ and SlewTime real , compares favorably to the worst-case bound of 2N .

Figure 1 .
Figure 1.Time dependent slew overheads at Keck Observatory.The left plot shows the motion of two targets in the alt/az frame across the accessible region at Keck Observatory (latitude = 19.8• N).The top target sits at δ = 43.8• , the bottom at δ = 11.7 • .As time progresses, the azimuthal separation dominates the slew and steadily increases to a maximum of over five minutes.Late in the night, the lower declination target crosses azimuth = 270 • , which corresponds to the assumed telescope cable wrap limit.At this time, a direct slew becomes available, and the slew time drops significantly.For any two arbitrary targets i and j, τ slew ijm approximates the pairwise slew time within the bounds of each time slot as dictated by M .The right plot shows the interpolated slew curve (grey) with the travel time tensor τ slew ijm (black) over-plotted for the full night with M = 3.In this extreme case, τ slew ijm often misrepresents the slew time by several minutes in the second and third time slot.

Figure 2 .
Figure2.Schematic of sub-optimal, optimal, and infeasible tours in the travelling telescope problem.The four polar plots on top show the sky above Keck Observatory (latitude = +19 deg) in alt/az coordinates.During a duration of four hours, targets A-E move across the sky as the Earth rotates.The shaded region shows an inaccessible region of alt/az (here, the Keck-I Naysmth platform).We have exaggerated slew times to clearly illustrate the differences between the three tours.In the top sub-optimal tour, we include a horizontal timeline for targets A-D.The windows of accessibility are shown as vertical ticks.For example, the earliest the telescope can complete observations of target A is labeled with t e A while t l A , denotes the latest time.Note that t e A depends on both the rise time and exposure duration, and t ℓ i is the either the set time or the end of the observing interval, whichever comes first.This scheduler uses a greedy approach; i.e., once the observation A of completes, the telescope immediately slews to the unobserved target with the shortest slew time.The tour is feasible, but the entire observing duration is required in this example.The middle optimal tour is scheduled with a global optimization.In our example, the slew between A and D is so long that it is advantageous for the telescope to wait until target B rises to observe it, before proceeding to D. This is the optimal tour to observe targets A-D.Finally the bottom plot shows the inclusion of a fifth target E to illustrate an infeasible tour.There is not enough time to observe all targets.

Figure 3 .
Figure 3.Comparison of a simple heuristic target ordering and a TTP optimized tour of N = 50 targets during a simulated half night.The former routine orders the targets by their increasing set times.Top row: the azimuthal coordinate of the simulated telescope in each case.Middle row: the elevation angle of the simulated telescope.Bottom row: time spent exposing, slewing, or idle.The targets have randomized time windows, but are all accessible for at least half of the scheduling duration.The worst case slew estimation from Equation11is 3.34 minutes per target.For the left (heuristic) column, the total slew time is 73 minutes, or an average of 91 seconds per slew, leaving 27 minutes idle.In the optimized (right) tour, the total time is a mere 13.8 minutes, or an average of 17.3 seconds per slew, resulting in 86 minutes idle (23 additional targets at the same pace).

Figure 4 .
Figure 4. Top row: average value of Runtime for TTP-Global when scheduling each duration type D, with varying values of M .In the quarter night case, TTP-Global can schedule all 25 observations for M < 10, but struggles to reach optimality at M = 10.In the half night case, only the M = 1 model achieves optimality for N = 50.For the full night, TTP-Global handles N = 100 far better than expected for the static case in comparison to N = 50.Bottom row: average relative improvement RelRedreal for each value of D across the model types.For all D and the respective maximum values of N , TTP-Global can reliably reduceslews by around a factor of 5.For small N , there is some benefit to using M > 1.For larger models, the risk of long slews in the M = 1 case has little impact on RelRedreal.
(Astropy Collaboration et al.   2013, 2018), numpy(Harris et al. 2020), pandas (pandas development team 2020), gurobi (Gurobi Optimization, LLC 2023), matplotlib(Hunter 2007) APPENDIX A. VARIABLES Table2lists the variables used in all preceding sections of this paper, along with their first usage: has the same structure as problem (C5), except that we force the violation variables ϵ slot,e z and ϵ slot,ℓ z to zero; thus, we no longer allow for violations of the slot bound constraints (C4h) and (C4i).We do still allow for violations of the rise and set time constraints (C4j) and (C4k).The objective function measures how much the rise and set time constraints are violated.Observe that if the objective value of problem (C6) is zero, then we have exactly verified that the minimum slew time MILP (C4) is feasible.

Figure 5 .
Figure 5. Runtime and RelRedreal for each duration type D, as found in Tables 1, 3, and 4. The first two rows summarize these statistics for the D = Quarter simulations, the next two for D = Half, and the bottom two for D = Full.The top row in each pair shows the average Runtime (and optimal subset, for TTP-Global) across all runs for D with varying M as a function of N .The bottom row in each pair shows the average RelRed real for the same values of D and M .

Table 2 .
Symbols UsedOne may require two exposures i and i ′ be scheduled back-to-back, such as a science observation and a calibration observation.Such linked observations may be specified via two additional constraints: