diff --git a/dorado/scheduling/mission.py b/dorado/scheduling/mission.py index 412e537..ea998d0 100644 --- a/dorado/scheduling/mission.py +++ b/dorado/scheduling/mission.py @@ -12,12 +12,14 @@ import astroplan from astropy import units as u +import numpy as np from . import data from .constraints import (BrightEarthLimbConstraint, EarthLimbConstraint, TrappedParticleFluxConstraint, get_field_of_regard) from .fov import FOV from .orbit import Orbit +from ._slew import slew_time __all__ = ('Mission', 'dorado', 'ultrasat', 'uvex') @@ -32,13 +34,34 @@ class Mission: """Container for mission configuration.""" constraints: Collection + """Observing constraint.""" + fov: FOV + """Field of view.""" + orbit: Orbit + """The orbit.""" + + min_overhead: u.Quantity + """Minimum overhead between observations (readout and settling time).""" + + max_angular_velocity: u.Quantity + """Maximum angular velocity for slews.""" + + max_angular_acceleration: u.Quantity + """Maximum angular acceleration for slews.""" def get_field_of_regard(self, *args, **kwargs): return get_field_of_regard(self.orbit, self.constraints, *args, **kwargs) + def overhead(self, coord1, coord2): + # FIXME: doesn't handle slews between different roll angles! + return np.maximum(self.min_overhead, + slew_time(coord1.separation(coord2), + self.max_angular_velocity, + self.max_angular_acceleration)) + dorado = Mission( constraints=( @@ -50,7 +73,10 @@ def get_field_of_regard(self, *args, **kwargs): astroplan.MoonSeparationConstraint(23 * u.deg), astroplan.GalacticLatitudeConstraint(10 * u.deg)), fov=FOV.from_rectangle(7.1 * u.deg), - orbit=_read_orbit('dorado-625km-sunsync.tle') + orbit=_read_orbit('dorado-625km-sunsync.tle'), + min_overhead=42 * u.s, + max_angular_velocity=0.872 * u.deg / u.s, + max_angular_acceleration=0.244 * u.deg / u.s**2 ) """Configuration for Dorado. @@ -66,6 +92,10 @@ def get_field_of_regard(self, *args, **kwargs): * The orbit one of four synthetic plausible orbits generated by Craig Markwardt. It is a nearly circular sun-synchronous orbit at an altitude of 625 km. + +* The maximum angular acceleration and angular velocity about its three + principal axes were given to us by the spacecraft vendor, but to be + conservative we are using the smallest values. """ @@ -76,7 +106,10 @@ def get_field_of_regard(self, *args, **kwargs): astroplan.MoonSeparationConstraint(23 * u.deg), astroplan.GalacticLatitudeConstraint(10 * u.deg)), fov=FOV.from_rectangle(14.1 * u.deg), - orbit=_read_orbit('goes17.tle') + orbit=_read_orbit('goes17.tle'), + min_overhead=42 * u.s, + max_angular_velocity=0.872 * u.deg / u.s, + max_angular_acceleration=0.244 * u.deg / u.s**2 ) """Configuration for ULTRASAT. @@ -94,6 +127,9 @@ def get_field_of_regard(self, *args, **kwargs): * ULTRASAT will be in a geosynchronous orbit. Here, we are using the orbital elements of GOES-17, a real geosynchronous weather satellite that is currently on orbit. + +* Maximum angular acceleration, maximum angular velocity, and overhead time are + assumed to be the same as for Dorado. """ @@ -104,7 +140,10 @@ def get_field_of_regard(self, *args, **kwargs): astroplan.MoonSeparationConstraint(23 * u.deg) ), fov=FOV.from_rectangle(3.3 * u.deg), - orbit=_read_orbit('vela1.tle') + orbit=_read_orbit('vela1.tle'), + min_overhead=42 * u.s, + max_angular_velocity=0.872 * u.deg / u.s, + max_angular_acceleration=0.244 * u.deg / u.s**2 ) """Configuration for UVEX. @@ -122,4 +161,7 @@ def get_field_of_regard(self, *args, **kwargs): with a perigee greater than 45,000 km. As a proxy, we use the actual two-line elements for VELA 1, a historically important gamma-ray satellite with a similar orbit. + +* Maximum angular acceleration, maximum angular velocity, and overhead time are + assumed to be the same as for Dorado. """ diff --git a/dorado/scheduling/scripts/main.py b/dorado/scheduling/scripts/main.py index 8641b29..a429940 100644 --- a/dorado/scheduling/scripts/main.py +++ b/dorado/scheduling/scripts/main.py @@ -85,39 +85,53 @@ def main(args=None): import os # import shlex import sys + import warnings from astropy_healpix import HEALPix from astropy.coordinates import ICRS from astropy.io import fits from astropy.time import Time from astropy.table import Table + from cplex.callbacks import LazyConstraintCallback from docplex.mp.model import Model + from docplex.mp.callbacks.cb_mixin import ConstraintCallbackMixin from ligo.skymap.io import read_sky_map from ligo.skymap.bayestar import rasterize from ligo.skymap.util import Stopwatch import numpy as np - from scipy.signal import convolve from tqdm import tqdm + from ..utils import nonzero_intervals + mission = getattr(_mission, args.mission) healpix = HEALPix(args.nside, order='nested', frame=ICRS()) log.info('reading sky map') # Read multi-order sky map and rasterize to working resolution - start_time = Time(fits.getval(args.skymap, 'DATE-OBS', ext=1)) + event_time = Time(fits.getval(args.skymap, 'DATE-OBS', ext=1)) skymap = read_sky_map(args.skymap, moc=True)['UNIQ', 'PROBDENSITY'] prob = rasterize(skymap, healpix.level)['PROB'] # Set up grids with u.add_enabled_equivalencies(equivalencies.orbital(mission.orbit)): - time_steps_per_exposure = int(np.round( - (args.exptime / args.time_step).to_value( - u.dimensionless_unscaled))) - times = start_time + args.delay + np.arange( - 0, args.duration.to_value(u.s), - args.time_step.to_value(u.s)) * u.s + nexp = int(((mission.min_overhead + args.duration) / + (mission.min_overhead + args.exptime)).to( + u.dimensionless_unscaled)) + if args.nexp is not None and args.nexp < nexp: + nexp = args.nexp + min_delay_s = (mission.min_overhead + args.exptime).to_value(u.s) + exptime_s = args.exptime.to_value(u.s) + duration_s = args.duration.to_value(u.s) + time_step_s = args.time_step.to_value(u.s) + start_time = event_time + args.delay + times = start_time + np.arange( + 0, duration_s, time_step_s) * u.s rolls = np.arange(0, 90, args.roll_step.to_value(u.deg)) * u.deg + if len(rolls) > 1: + warnings.warn( + 'Slew constraints for varying roll angles are not yet implemented') + if args.skygrid_file is not None: centers = Table.read(args.skygrid_file, format='ascii.ecsv')['center'] else: @@ -125,10 +139,7 @@ def main(args=None): args.skygrid_step) log.info('evaluating field of regard') - not_regard = convolve( - ~mission.get_field_of_regard(centers, times, jobs=args.jobs), - np.ones(time_steps_per_exposure)[:, np.newaxis], - mode='valid', method='direct') + regard = mission.get_field_of_regard(centers, times, jobs=args.jobs) log.info('generating model') m = Model() @@ -137,44 +148,41 @@ def main(args=None): if args.jobs is not None: m.context.cplex_parameters.threads = args.jobs - log.info('adding variable: observing schedule') - shape = (len(centers), len(rolls), not_regard.shape[0]) - schedule = np.reshape(m.binary_var_list(np.prod(shape)), shape) + log.info('variable: start time for each observation') + obs_start_time = np.asarray( + m.continuous_var_list(nexp, lb=0, ub=duration_s - exptime_s)) - log.info('adding variable: whether a given pixel is observed') - pixel_observed = np.asarray(m.binary_var_list(healpix.npix)) + log.info('variable: field selection for each observation') + shape = (nexp, len(centers), len(rolls)) + obs_field = np.reshape(m.binary_var_list(np.prod(shape)), shape) - log.info('adding variable: whether a given field is used') - field_used = np.reshape(m.binary_var_list(np.prod(shape[:2])), shape[:2]) + log.info('variable: whether a given observation is used') + obs_used = np.asarray(m.binary_var_list(nexp)) - log.info('adding variable: whether a given time step is used') - time_used = np.asarray(m.binary_var_list(shape[2])) + log.info('variable: whether a given field is used') + field_used = np.reshape(m.binary_var_list(np.prod(shape[1:])), shape[1:]) - if args.nexp is not None: - log.info('adding constraint: number of exposures') - m.add_constraint_(m.sum(time_used) <= args.nexp) + log.info('variable: whether a given pixel is observed') + pix_used = np.asarray(m.binary_var_list(healpix.npix)) - log.info('adding constraint: only observe one field at a time') - m.add_constraints_( - m.sum(schedule[..., i].ravel()) <= 1 for i in tqdm(range(shape[2])) - ) + log.info('constraint: at most one field is used for each observation') + for i in tqdm(range(nexp)): + m.add_sos1(obs_field[i].ravel()) + + log.info('constraint: consecutive observations are used') + m.add_(obs_used[i] >= obs_used[i + 1] for i in tqdm(range(nexp - 1))) + + log.info('constraint: an observation is used if any field is chosen') m.add_equivalences( - time_used, - [m.sum(schedule[..., i].ravel()) >= 1 for i in tqdm(range(shape[2]))] - ) - m.add_constraints_( - m.sum(time_used[i:i+time_steps_per_exposure]) <= 1 - for i in tqdm(range(schedule.shape[-1])) - ) + obs_used, + (m.logical_or(*obs_field[i].ravel()) >= 1 for i in tqdm(range(nexp)))) - log.info('adding constraint: a pixel is observed if it is in any field') - m.add_constraints_( - m.sum(lhs) >= rhs - for lhs, rhs in zip( - tqdm(schedule.reshape(field_used.size, -1)), - field_used.ravel() - ) - ) + log.info('constraint: a field is used if it is chosen for an observation') + m.add_( + m.sum(obs_field[:, j, k]) >= field_used[j, k] + for j in tqdm(range(len(centers))) for k in range(len(rolls))) + + log.info('constraint: a pixel is used if it is in any used fields') indices = [[] for _ in range(healpix.npix)] with tqdm(total=len(centers) * len(rolls)) as progress: for i, grid_i in enumerate( @@ -183,17 +191,98 @@ def main(args=None): for k in grid_ij: indices[k].append((i, j)) progress.update() - m.add_constraints_( + m.add_( m.sum(field_used[lhs_index] for lhs_index in lhs_indices) >= rhs - for lhs_indices, rhs in zip(tqdm(indices), pixel_observed) - ) - - log.info('adding constraint: field of regard') - i, j = np.nonzero(not_regard) - m.add_constraint_(m.sum(schedule[j, :, i].ravel()) <= 0) + for lhs_indices, rhs in zip(tqdm(indices), pix_used)) + + log.info('constraint: observations do not overlap') + m.add_indicator_constraints_( + obs_used[i + 1] >> + (obs_start_time[i + 1] - obs_start_time[i] >= min_delay_s) + for i in range(nexp - 1)) + + log.info('constraint: field of regard') + for j in tqdm(range(len(centers))): + for k in range(len(rolls)): + # FIXME: not roll dependent yet + intervals = nonzero_intervals(regard[:, j]) * time_step_s + + if len(intervals) == 0: + # The field is always outside the field of regard, + # so disallow it entirely. + m.add_(m.sum(obs_field[:, j, k]) <= 0) + elif len(intervals) == 1: + # The field is within the field of regard during a single + # contigous interval, so require the observaiton to be within + # that interval. + (interval_start, interval_end), = intervals + m.add_indicator_constraints_( + obs_field[i, j, k] >> + (obs_start_time[i] >= interval_start) + for i in range(nexp)) + m.add_indicator_constraints_( + obs_field[i, j, k] >> + (obs_start_time[i] <= interval_end - exptime_s) + for i in range(nexp)) + else: + # The field is within the field of regard during two or more + # disjoint intervals, so introduce additional decision + # variables to decide which interval. + interval_start, interval_end = intervals.T + shape = (nexp, len(intervals)) + interval_choice = np.reshape( + m.binary_var_list(np.prod(shape)), shape) + m.add_indicator_constraints_( + obs_field[i, j, k] >> (m.sum(interval_choice[i]) >= 1) + for i in range(nexp)) + m.add_indicator_constraints_( + interval_choice[i, i1] >> + (obs_start_time[i] >= interval_start[i1]) + for i in range(nexp) for i1 in range(len(intervals))) + m.add_indicator_constraints_( + interval_choice[i, i1] >> + (obs_start_time[i] <= interval_end[i1] - exptime_s) + for i in range(nexp) for i1 in range(len(intervals))) + + class SlewConstraintCallback(ConstraintCallbackMixin, + LazyConstraintCallback): + + def __init__(self, env): + LazyConstraintCallback.__init__(self, env) + ConstraintCallbackMixin.__init__(self) + + def __call__(self): + # Reconstruct partial solution + sol = self.make_solution_from_watched() + + # Determine which fields are selected + i, j, k = np.nonzero(np.reshape(sol.get_values(obs_field.ravel()), + obs_field.shape)) + + # Calculate overhead + exptime between each pair of fields + coords = centers[j] + dt_array = mission.overhead(coords[1:], coords[:-1]).to_value(u.s) + dt_array += exptime_s + + lhs_array = obs_start_time[i] + rhs_array = obs_field[i, j, k] + # This is a big-M formulation of the indicator constraint: + # (rhs1 & rhs0) >> lhs1 - lhs0 >= duration_s + cons = [ + lhs1 - lhs0 - dt >= duration_s * (rhs1 + rhs0 - 2) + for lhs0, lhs1, rhs0, rhs1, dt + in zip(lhs_array[:-1], lhs_array[1:], + rhs_array[:-1], rhs_array[1:], dt_array)] + + for _, lhs, sense, rhs in self.get_cpx_unsatisfied_cts(cons, sol): + self.add(lhs, sense, rhs) + + cb = m.register_callback(SlewConstraintCallback) + cb.register_watched_vars(obs_field.ravel()) + cb.register_watched_vars(obs_start_time) log.info('adding objective') - m.maximize(m.scal_prod(pixel_observed, prob)) + m.maximize(m.scal_prod(pix_used, prob)) log.info('solving') stopwatch = Stopwatch() @@ -203,24 +292,27 @@ def main(args=None): log.info('extracting results') if solution is None: - schedule_flags = np.zeros(schedule.shape, dtype=bool) + obs_field_value = np.zeros(obs_field.shape, dtype=bool) + obs_start_time_value = np.zeros(obs_start_time.shape) objective_value = 0.0 else: - schedule_flags = np.asarray( - solution.get_values(schedule.ravel()), dtype=bool + obs_field_value = np.asarray( + solution.get_values(obs_field.ravel()), dtype=bool ).reshape( - schedule.shape + obs_field.shape ) + obs_start_time_value = np.asarray(solution.get_values(obs_start_time)) objective_value = m.objective_value + obs_start_time_value = start_time + obs_start_time_value * u.s - ipix, iroll, itime = np.nonzero(schedule_flags) + i, j, k = np.nonzero(obs_field_value) result = Table( data={ - 'time': times[itime], - 'exptime': np.repeat(args.exptime, len(times[itime])), - 'location': mission.orbit(times).earth_location[itime], - 'center': centers[ipix], - 'roll': rolls[iroll] + 'time': obs_start_time_value[i], + 'exptime': np.repeat(args.exptime, len(obs_start_time_value[i])), + 'location': mission.orbit(obs_start_time_value).earth_location[i], + 'center': centers[j], + 'roll': rolls[k] }, descriptions={ 'time': 'Start time of observation',