from __future__ import division, print_function
import os, shutil, subprocess, sys, time, glob
import pandas as pd
import openpyxl
import pyomo.environ as po
from pyomo.opt import SolverStatus, TerminationCondition

# choose setpoints to optimize HVAC operation for one month, based on prices in
# pjmLMP and weather and regression coefficients in Excel workbook

# will need to save settings in temperature_setpoints.tsv (one year of hourly
# setpoints for 5 zones, with header, tab-separated), then execute
# runenergyplus RefBldgLargeOfficeNew2004_Chicago_August_Scheduled.idf USA_IL_Chicago-OHare.Intl.AP.725300_TMY3.epw

# Then export 5-minute values for all zones and use the
# same ipynb script to read those along with the MLEP results and show the results.

# TODO: setup a schedule with random up-down jumps in setpoints, then use
# that as the basis for the regression and see if you get more realistic
# coefficients (gain equation could just be wrong too). Right now, the model
# thinks it has to run at full power all day to keep the place cool. Or it does
# too much cooling at night and wastes money (with poor comfort too, because
# it chooses high setpoints because it thinks low ones are unreachable).

# constants for all methods
t_ideal = 22.5
n_floors = 11  # number of repetitions of each floor's population
energy_plus_dir = 'RefLargeOfficeBuilding'
timestep_building = 'RefBldgLargeOfficeNew2004_Chicago_timestep.idf'
hourly_building =  'RefBldgLargeOfficeNew2004_Chicago_hourly.idf'

# note: for the pso runs, Energy Plus was setup to use daylight saving time;
# in this condition, all inputs and outputs are still done in standard time,
# (and I think the EMS controls we added also work in standard time),
# but built-in occupancy and equipment schedules get moved one hour earlier.
# However, PSO was optimized as if the occupancy schedules matched the input-
# output times (which made the PSO program assume the building was occupied
# one hour later than Energy Plus would have reported). For consistency, we
# use the same occupancy schedule here that was used for PSO, rather than
# retrieving occupancy data from Energy Plus.
# data from Matlab PSO program:
weekdayEndHours = [0, 8, 11, 12 ,13 ,14 ,18, 19, 24]
weekdayOccupancy = [0, 0, 1,  .8, .4, .8, 1,  .5, 0]
weekendEndHours = [0, 8, 9,   11  ,13 ,16 ,18, 19, 24]
weekendOccupancy = [0, 0, 0.2, 0.3, 0, .3, .2,  .1, 0]
occupancy = pd.Series(index=pd.date_range('2007-08-01 00:00:00', '2007-08-31 23:55:00', freq='5min'))
dt = occupancy.index
for i in range(1, len(weekdayEndHours)):
    occupancy[
        (dt.dayofweek < 5)
        & (dt.hour >= weekdayEndHours[i-1])
        & (dt.hour < weekdayEndHours[i])
    ] = weekdayOccupancy[i]
    occupancy[
        (dt.dayofweek >= 5)
        & (dt.hour >= weekendEndHours[i-1])
        & (dt.hour < weekendEndHours[i])
    ] = weekendOccupancy[i]

# TO USE:
# import this module in Python (or open it with Atom editor using Hydrogen package)
# 1. run calibrate_leow and then follow the instructions there to configure the
# leow model.
# 2. run all the run_* functions to produce detailed and summary outputs (note
# that run_ccpso_results is only used for double-checking; we report results
# from the original ccpso run, placed in energy_stats.csv by aggregate_energy_data7.ipynb
# 3. open the summary*.csv files in Excel and paste the data into
# "results summary and graphs original outputs.xlsx". 

def calibrate_leow():
    setpoints = get_night_setback_setpoints(22.5)
    # setpoints = get_random_setpoints()
    runenergyplus(hourly_building, setpoints, 'calibration_data.csv')
    # copy calibration_data.csv into Excel linear model file
    # (could be streamlined to load the results, run a linear regression
    # and save the regression coefficients for later reading by get_leow_setpoints())

def run_leow_optimization():
    results = []
    for penalty in 0.01 * 10**pd.np.linspace(0, -1.699, 16):
        setpoints = get_leow_setpoints(float(penalty))
        c, d = runenergyplus(timestep_building, setpoints, 'leow_opt_{}.csv'.format(penalty))
        results.append((str(penalty), c, d))
    results_df = pd.DataFrame.from_records(results, columns=['scenario', 'cost', 'discomfort'])
    results_df.to_csv(os.path.join('results', 'leow_summary.csv'), index=False)

def run_night_setback():
    results = []
    for t_day in pd.np.arange(22.5, 24.5, 0.125):
        setpoints = get_night_setback_setpoints(float(t_day))
        c, d = runenergyplus(timestep_building, setpoints, 'night_setback_{}.csv'.format(t_day))
        results.append((str(t_day), c, d))
    results_df = pd.DataFrame.from_records(results, columns=['scenario', 'cost', 'discomfort'])
    results_df.to_csv(os.path.join('results', 'night_setback_summary.csv'), index=False)

if False:
    t_day = 22.5
    setpoints = get_night_setback_setpoints(float(t_day))
    runenergyplus(timestep_building, setpoints, 'night_setback_on_sun_{}.csv'.format(t_day))

def run_transactive_power():
    # values from table 3.4 on p. 3.15 of PNNL report (converted from deg F to deg C):
    # https://www.pnnl.gov/main/publications/external/technical_reports/PNNL-17167.pdf)
    settings = [
        (k, dt_min, dt_max)
        for dt_min in [0, -3/1.8]
        for dt_max in [10/1.8, 5/1.8]
        for k in [1, 2, 3]
    ]
    results = []
    for i, (k, delta_t_min, delta_t_max) in enumerate(settings):
        setpoints = get_transactive_power_setpoints(k, delta_t_min, delta_t_max)
        c, d = runenergyplus(timestep_building, setpoints, 'transactive_power_{}.csv'.format(i))
        results.append((str(i), c, d))
    results_df = pd.DataFrame.from_records(results, columns=['scenario', 'cost', 'discomfort'])
    results_df.to_csv(os.path.join('results', 'transactive_power_summary.csv'), index=False)


def run_ccpso_results():
    results = []
    for path in glob.glob(os.path.join('..', 'final', 'Window_Shifting_Logs', '*')):
        tag = '_'.join(os.path.split(path)[-1].split('_')[1:])
        setpoints = get_ccpso_setpoints(path)
        if tag.startswith('Corbin'):
            file = '{}.csv'.format(tag)
        else:
            file = 'ccpso_{}.csv'.format(tag)
        c, d = runenergyplus(timestep_building, setpoints, file)
        results.append((tag, c, d))
    results_df = pd.DataFrame.from_records(results, columns=['scenario', 'cost', 'discomfort'])
    results_df.to_csv(os.path.join('results', 'ccpso_summary.csv'), index=False)


def get_power_consumption_statistics():
    scenarios = [
        'night_setback_24.0.csv', 'leow_opt_0.000567631733412.csv',
        'transactive_power_2.csv', 'Corbin_WithTerm_10e-3.csv', 'ccpso_1.9e-3.csv'
    ]
    prices = get_prices()
    results = pd.DataFrame()
    daily_energy = pd.DataFrame()
    for scenario in scenarios:
        tag = scenario[:-4].replace('_', ' ')
        df = pd.read_csv(os.path.join('results', scenario))
        dt = pd.Index(get_energyplus_datetime_index(df))
        hour = dt - pd.to_timedelta(dt.minute, 'm')
        df['price'] = prices[hour].values
        # power in kW
        df['power'] = df['Whole Building:Facility Total HVAC Electric Demand Power [W](TimeStep)']*0.001
        low_price = df['price'].quantile(0.05)
        high_price = df['price'].quantile(0.95)
        # get daily power use in MWh
        daily_energy[tag] = df.groupby(dt.date)['power'].mean() * 24 * 0.001

        stats = pd.DataFrame(
            dict(
                cost=(df['power']*df['price']/1000/12).sum(),
                avg_load_kw=df['power'].mean(),
                avg_load_top=df.loc[df['price']>=high_price, 'power'].mean(),
                avg_load_bottom=df.loc[df['price']<=low_price, 'power'].mean(),
                avg_price=(df['power']*df['price']).sum()/df['power'].sum()
            ),
            index=[tag]
        )
        results = results.append(stats)

    results.to_csv(os.path.join('results', 'cost_statistics.csv'))
    daily_energy.to_csv(os.path.join('results', 'daily_energy.csv'))

def get_energyplus_datetime_index(df):
    """ Retrieve Date/Time column from dataframe df and return as a datetime index (assumed to be in 2007). """
    dt = (
        pd.to_datetime(df['Date/Time'].str[1:6] + '/2007')
        + pd.to_timedelta(df['Date/Time'].str[7:])
    )
    # step back by one timestep, since Energy Plus uses end-of-step indexing
    dt -= (dt[1] - dt[0])
    return dt


def runenergyplus(building, setpoints, output_file):
    # (building, setpoints, output_file) = (timestep_building, setpoints, 'leow_opt_{}.csv'.format(penalty))
    # make a setpoint schedule file
    save_setpoints(
        setpoints,
        os.path.join(energy_plus_dir, 'temperature_setpoints.csv')
    )
    # call runenergyplus with specified building in energy_plus_directory
    # (will read the saved setpoints)
    print("Running EnergyPlus...")
    start = time.time()
    status = subprocess.call(
        ['runenergyplus', building, 'USA_IL_Chicago-OHare.Intl.AP.725300_TMY3.epw'],
        cwd=energy_plus_dir
    )
    print("EnergyPlus finished in {} seconds.".format(time.time()-start))
    if status != 0:
        raise RuntimeError(
            'EnergyPlus finished with error code {}. See {} for details.'
            .format(
                status,
                os.path.join(energy_plus_dir, 'Output', building.replace('.idf', '.err'))
            )
        )
    # save results for later reference
    if not os.path.exists('results'):
        os.makedirs('results')
    out_file = os.path.join('results', output_file)
    shutil.copyfile(
        os.path.join(energy_plus_dir, 'Output', building.replace('.idf', '.csv')),
        out_file
    )

    print("Results stored in {}.".format(out_file))
    if building == timestep_building:
        return get_cost_and_discomfort(out_file)
    else:
        return None

def get_cost_and_discomfort(file):
    # file='results/ccpso_1.7e-3.csv'
    # file='results/ccpso_3.2e-3.csv'
    zone_pops = data_frame_from_xlsx('Linear Building Model.xlsx', 'zone_populations').T.set_index(0).T.astype(float)
    zone_pops.columns = zone_pops.columns.str.replace('Hourly', 'TimeStep')
    prices = get_prices()
    results = pd.read_csv(file)
    # create index, including missing year and converting from 00:00:05-24:00:00 to 00:00:00-23:55:00
    results.index = get_energyplus_datetime_index(results)

    zone_temps = results.loc[:, zone_pops.columns].values
    population = occupancy.to_frame().values * zone_pops.values * n_floors

    # calculate total discomfort, including C-F conversion and adjustment for
    # number of timesteps per hour (12)
    total_discomfort = ((((zone_temps-t_ideal)*1.8)**2) * population).sum() / 12

    power = results['Whole Building:Facility Total HVAC Electric Demand Power [W](TimeStep)']
    power.index = power.index - pd.to_timedelta(power.index.minute, 'm')
    total_cost = (power * prices).dropna().sum() / 1e6 / 12

    return (total_cost, total_discomfort)

def save_setpoints(setpoints, file):
    """
    Save setpoints to the specified file; setpoints should be a pandas series
    with a datetime index spanning dates in the year 2007; t_ideal will be used
    for any missing setpoints.
    """
    settings = pd.DataFrame(index=pd.date_range('2007', periods=8760, freq='h'))
    for i in range(1, 6):
        settings['Tset{}'.format(i)] = setpoints
    # settings.loc['2007-08-15':'2007-08-18', :]
    settings = settings.fillna(22.5)
    # Energy Plus uses the first entry as the last hour of the same day (treats
    # it as ending at midnight at the end of the day, even though it's sort of
    # in position to end at midnight at the start of the day), then goes through
    # the other 23 hours. So we have to re-sort the hours to get them used in
    # the right order.
    settings['date'] = settings.index.date
    settings['hour'] = settings.index.hour
    settings['hour'] = (settings['hour'] + 1).mod(24)
    # or settings['hour'] = (settings.index.to_series().dt.hour+1).mod(24)

    settings = settings.sort_values(['date', 'hour'])
    settings.to_csv(file, index=True)


def get_ccpso_setpoints(path):
    """
    retrieve setpoints previously stored in
    path/CSV_Data-Planned_{31-Aug-2012_DOY243}/Day70.csv
    and aggregate into a standard setpoints series.
    """
    ts = pd.Series([0] * 24*31)
    ts.index = pd.date_range('2007-08-01', periods=len(ts), freq='h')

    for date in pd.date_range('2007-08-01', periods=31, freq='d'):
        subdir = 'CSV_Data-Planned_' +  date.strftime('%d-%b-%Y_DOY%j').replace('2007', '2012')
        day = pd.read_csv(os.path.join(path, subdir, 'Day70.csv'))
        day.index = pd.date_range(date, periods=24*12, freq='300s')
        day = day[day.index.minute == 0]
        ts[day.index] = day['Tset_Zone1 (C)']

    return ts

def aggregate_original_results():
    # save values for cross-checking with rerun results
    path = '../final/Window_Shifting_Logs/HorizonLogs_1.7e-3'
    df = pd.DataFrame()  # index = pd.date_range('2007-08-01', periods=len(ts), freq='h'))
    for date in pd.date_range('2007-08-01', periods=31, freq='d'):
        subdir = 'CSV_Data-Planned_' +  date.strftime('%d-%b-%Y_DOY%j').replace('2007', '2012')
        day = pd.read_csv(os.path.join(path, subdir, 'Day70.csv'))
        day.index = pd.date_range(date, periods=24*12, freq='300s')
        df = df.append(day)

    # save values for cross-checking with rerun results
    df.to_csv('results/original_values_1.7e-3.csv')

def create_rerun_setpoints():
    path = '../final/Window_Shifting_Logs/HorizonLogs_1.7e-3'
    df = pd.DataFrame()  # index = pd.date_range('2007-08-01', periods=len(ts), freq='h'))
    for date in pd.date_range('2007-08-01', periods=31, freq='d'):
        subdir = 'CSV_Data-Planned_' +  date.strftime('%d-%b-%Y_DOY%j').replace('2007', '2012')
        day = pd.read_csv(os.path.join(path, subdir, 'Day70.csv'))
        day.index = pd.date_range(date, periods=24*12, freq='300s')
        day = day[day.index.minute == 0]
        df = df.append(day)

    # make a full year and re-sort to use the right setpoints at the right times
    df = df.reindex(pd.date_range('2007', periods=8760, freq='h'))
    df = df.fillna(0.0)
    # Energy Plus uses the first entry as the last hour of the same day (treats
    # it as ending at midnight at the end of the day, even though it's sort of
    # in position to end at midnight at the start of the day), then goes through
    # the other 23 hours. So we have to re-sort the hours to get them used in
    # the right order.
    df['date'] = df.index.date
    df['hour'] = df.index.hour
    df['hour'] = (df['hour'] + 1).mod(24)

    df = df.sort_values(['date', 'hour'])
    df.to_csv('RefLargeOfficeBuilding/rerun_setpoints.csv', index=True)


# leow optimization will need to retrieve hourly outdoor temperatures from weather file...
def get_leow_setpoints(penalty=0.038, use_setback=False):
    prices = get_prices()
    # use string indexes to work more easily with pyomo
    # prices.index = prices.index.to_series().astype(str)

    # get weather from regression file
    weather = data_frame_from_xlsx('Linear Building Model.xlsx', 'simulation_data').T.set_index(0).T
    # In EnergyPlus outputs, all data are shown as averages for the preceding hour,
    # so 1:00 means the hour from 0:00 to 1:00.
    # Add missing year and convert from 01:00:00-24:00:00 to 00:00:00-23:00:00
    weather['date_time'] = get_energyplus_datetime_index(weather)
    weather = weather.drop('Date/Time', axis=1)
    weather = weather.set_index('date_time')
    weather = weather.astype(float)

    # regression coefficients
    coeffs = data_frame_from_xlsx('Linear Building Model.xlsx', 'regression_coefficients').T.set_index(0).iloc[:, 0]

    # zone populations
    zone_pops = data_frame_from_xlsx('Linear Building Model.xlsx', 'zone_populations').T.set_index(0).iloc[:, 0].astype(float)

    # setup optimization model
    m = po.ConcreteModel()
    m.HOURS = po.Set(initialize=weather.index, ordered=True)
    m.Energy = po.Var(m.HOURS, within=po.NonNegativeReals)
    m.Temperature = po.Var(m.HOURS)
    m.occupancy = po.Param(m.HOURS, initialize=lambda m, h: occupancy[h])
    m.outdoor_t0 = po.Param(m.HOURS, initialize=lambda m, h: weather.loc[h, 'Outdoor T0'])
    m.electricity_price = po.Param(m.HOURS, initialize=lambda m, h: prices[h])
    m.max_people = po.Param(initialize=zone_pops.sum() * 11) # account for having 11 floors
    m.discomfort_penalty = penalty  # $ per person per degree C ** 2 per hour
    m.t_ideal = t_ideal
    m.Cost = po.Objective(
        expr=sum(
            m.Energy[h] * m.electricity_price[h]
            + m.discomfort_penalty * m.max_people * m.occupancy[h] * (m.Temperature[h] - m.t_ideal) ** 2
            for h in m.HOURS
        ),
        sense=po.minimize
    )
    m.System_Size = po.Constraint(m.HOURS, rule=lambda m, h: m.Energy[h] <= 1.0)
    m.Building_Model = po.Constraint(m.HOURS, rule=lambda m, h:
        m.Temperature[h]
        ==
        coeffs['Prev Hour T'] * m.Temperature[m.HOURS.prevw(h)]
        + coeffs['Outdoor T0'] * m.outdoor_t0[h]
        + coeffs['Energy in (MWh)'] * m.Energy[h]
    )
    if use_setback:
        m.Night_Setback = po.Constraint(m.HOURS, rule=lambda m, h:
            po.Constraint.Skip
            if m.occupancy[h] > 0 or m.occupancy[m.HOURS.nextw(h, 3)] > 0
            else (m.Energy[h] == 0.0)
        )

    solve_pyomo_model(m)

    cols = [
        'Temperature', 'Energy', 'outdoor_t0', 'electricity_price', 'occupancy'
    ]
    results = pd.DataFrame(
        [
            tuple(po.value(getattr(m, c)[h]) for c in cols)
            for h in m.HOURS
        ],
        columns=cols
    )
    results.index = weather.index

    if use_setback:
        # force Tset high to get coasting during night setback periods
        for h in results.index:
            if m.occupancy[h] > 0 or m.occupancy[m.HOURS.nextw(h, 3)] > 0:
                pass
            else:
                results.loc[h, 'Temperature'] = 27

    # %matplotlib inline
    # results.loc['2007-08-01':'2007-08-31', 'Energy'].plot()
    # results.loc['2007-08-01':'2007-08-31', 'Temperature'].plot()
    # results[results['Energy']<0.01]
    # results.loc['2007-08-07':'2007-08-08']

    return results['Temperature']


def get_prices():
    """
    retrieve hourly electricity prices for the study period
    """
    # dates are pretty weird in this file -- 2013 for the first half of the year, then 2007
    # Note: this uses central standard time at start and end of year and
    # daylight saving time in the middle; it is also missing a record for
    # 3/10/13 02:00 (start of daylight saving) and has a non-included extra record
    # at the end of daylight saving time. These shouldn't cause problem, since we
    # just use a strip in the summer which is all in CDT.
    prices = pd.read_excel('pjmLMP.xlsx')
    prices = prices.set_index('PUBLISHDATE').loc[:, 'H1':'H24']
    prices.index.name = 'date'
    prices.columns = list(range(24)) # switch to 0-based hours
    prices.columns.name = 'hour'
    prices = prices.stack().reset_index(name='price')
    prices['date_time'] = prices['date'] + pd.to_timedelta(prices['hour'], 'H')
    # prices['idx'] = prices['date_time'].astype(str)
    prices = prices.set_index('date_time')['price']
    return prices


def get_night_setback_setpoints(t_day=t_ideal):
    """
    Generate setpoints that follow a night-setback schedule for the
    RefBldgLargeOfficeNew2004 building (we have set this to have some occupancy
    at 8am to 7pm on all days)
    (could be updated to implement the approach we used for the HVAC paper instead)
    TODO: update to account for Sunday occupancy, i.e., follow our occupancy numbers,
    not the standard ones for this building
    """
    ts = pd.Series([26.7] * 8760)  # night schedule
    ts.index = pd.date_range('2007', periods=len(ts), freq='h')
    dt = ts.index

    # follow occupancy patterns (could give 1 hour head start, but then
    # this would differ from the extreme-cooling strategy used to evaluate PSO
    # termination period)
    ts[(dt.hour >= 8) & (dt.hour < 19)] = t_day
    # BLDG_OCC_SCH,            !- Name
    # For: Weekdays,           !- Field 9
    # Until: 8:00, 0.00,       !- Field 11
    # Until: 11:00, 1.00,      !- Field 13
    # Until: 12:00, 0.80,      !- Field 15
    # Until: 13:00, 0.40,      !- Field 17
    # Until: 14:00, 0.80,      !- Field 19
    # Until: 18:00, 1.00,      !- Field 21
    # Until: 19:00, 0.50,      !- Field 23
    # Until: 21:00, 0.00,      !- Field 25
    # Until: 24:00, 0.00,      !- Field 27
    # For: Weekends,           !- Field 28
    # Until: 8:00, 0.00,       !- Field 11
    # Until: 9:00, 0.10,      !- Field 13
    # Until: 11:00, 0.30,      !- Field 15
    # Until: 13:00, 0.00,      !- Field 17
    # Until: 16:00, 0.30,
    # Until: 18:00, 0.20,      !- Field 19
    # Until: 19:00, 0.10,      !- Field 21
    # Until: 24:00, 0.00,      !- Field 27

    # follow original cooling setpoint
    # ts[(dt.weekday <= 4) & (dt.hour >= 6) & (dt.hour < 22)] = t_day
    # ts[(dt.weekday == 5) & (dt.hour >= 6) & (dt.hour < 18)] = t_day
    # CLGSETP_SCH,             !- Name
    # Temperature,             !- Schedule Type Limits Name
    # Through: 12/31,          !- Field 1
    # For: Weekdays SummerDesignDay,  !- Field 2
    # Until: 06:00, 26.7,      !- Field 4
    # Until: 22:00, 24.0,      !- Field 6
    # Until: 24:00, 26.7,      !- Field 8
    # For: Saturday,           !- Field 9
    # Until: 06:00, 26.7,      !- Field 11
    # Until: 18:00, 24.0,      !- Field 13
    # Until: 24:00, 26.7,      !- Field 15
    # For WinterDesignDay,     !- Field 16
    # Until: 24:00, 26.7,      !- Field 18
    # For: AllOtherDays,       !- Field 19
    # Until: 24:00, 26.7;      !- Field 21
    return ts


def get_random_setpoints():
    import random
    ts = pd.Series([0.0] * 8760)
    t = t_ideal
    upper_lim = t_ideal + 2.5
    lower_lim = t_ideal - 2.5
    for i in range(len(ts)):
        adj = random.normalvariate(0, 1)
        if abs(adj) > 4:
            adj = 0
        t += adj
        if t < lower_lim:
            t = 2 * lower_lim - t
        if t > upper_lim:
            t = 2 * upper_lim - t
        ts[i] = t
    ts.index = pd.date_range('2007', periods=len(ts), freq='h')
    # %matplotlib inline
    # ts.plot()
    return ts

    settings = pd.DataFrame()
    for i in range(1, 6):
        settings['Tset{}'.format(i)] = ts
    %matplotlib inline
    settings.plot()

    # Energy Plus uses the first entry as the last hour of the same day (treats
    # it as ending at midnight at the end of the day, even though it's sort of
    # in position to end at midnight at the start of the day), then goes through
    # the other 23 hours. So we have to re-sort the hours to get them used in
    # the right order.
    settings['date'] = settings.index.date
    settings['hour'] = (settings.index.to_series().dt.hour+1).mod(24)
    settings = settings.sort_values(['date', 'hour'])
    settings.to_csv('temperature_setpoints_random.csv', index=True)



def solve_pyomo_model(m):
    solver = po.SolverFactory('cplexamp', solver_io='nl')

    print("Solving pyomo model...")
    results = solver.solve(m, tee=True, options_string='display=1')
    m.solutions.load_from(results)

    if results.solver.termination_condition == TerminationCondition.optimal:
        print("Solved model.")
        print("  Solver Status: " + str(results.solver.status))
        print("  Termination Condition: " + str(results.solver.termination_condition))
        return
    elif (results.solver.termination_condition == TerminationCondition.infeasible):
        if hasattr(model, "iis"):
            print("Model was infeasible; irreducibly inconsistent set (IIS) returned by solver:")
            print("\n".join(c.name for c in model.iis))
        else:
            print("Model was infeasible; if the solver can generate an irreducibly inconsistent set (IIS),")
            print("more information may be available by setting the appropriate flags in the ")
            print('solver_options_string and calling this script with "--suffixes iis".')
        raise RuntimeError("Infeasible model")
    else:
        print("Solver terminated abnormally.")
        print("  Solver Status: ", results.solver.status)
        print("  Termination Condition: ", results.solver.termination_condition)
        raise RuntimeError("Solver failed to find an optimal solution.")

def get_transactive_power_setpoints(k, delta_t_min, delta_t_max):
    """
    interpolate between highest
    and lowest comfortable temperature during occupied times, based on price;
    apply this to last two hours before occupancy if precooling is specified
    """
    # precool = True; k = 1
    max_occ_temp = t_ideal + delta_t_max
    min_occ_temp = t_ideal - delta_t_min

    setpoints = get_night_setback_setpoints(t_ideal)
    if delta_t_min > 0:  # precooling allowed
        setpoints[(setpoints.index.hour >= 6) & (setpoints.index.hour <= 7)] = t_ideal
    on_hours = (setpoints == t_ideal)

    # interpolate between ideal price and limits based on the prices during on_hours
    prices = get_prices()

    df = setpoints.to_frame(name='Tset')
    df = df.join(prices)  # omits first half of the year, which has 2013 prices

    # We use mean stddev of price during on_hours, not during last 24 hours,
    # so we remove prices for out-of-service hours
    df.loc[df['Tset']!=t_ideal, 'price'] = pd.np.nan
    df['mean'] = df.groupby(df.index.date)['price'].transform(lambda x: x.mean())
    df['std'] = df.groupby(df.index.date)['price'].transform(lambda x: x.std())

    # get active bound on temperature
    df['t_limit'] = max_occ_temp
    df.loc[df['price']<df['mean'], 't_limit'] = min_occ_temp

    # set temperature during occupied and pre-cooling periods, based on equ 2.2, p. 2.7 of
    # https://www.pnnl.gov/main/publications/external/technical_reports/PNNL-17167.pdf
    df['t_target'] = (
        t_ideal
        + (df['price'] - df['mean'])
        * (df['t_limit']-t_ideal).abs()
        / (k * df['std'])
    )
    df.loc[on_hours, 'Tset'] = df.loc[on_hours, 't_target']
    return df['Tset']


def data_frame_from_xlsx(xlsx_file, range_name):
    """ Get a single rectangular region from the specified file.
    range_name can be a standard Excel reference ('Sheet1!A2:B7') or
    refer to a named region ('my_cells')."""
    wb = openpyxl.load_workbook(xlsx_file, data_only=True, read_only=True)
    if '!' in range_name:
        # passed a worksheet!cell reference
        ws_name, reg = range_name.split('!')
        if ws_name.startswith("'") and ws_name.endswith("'"):
            # optionally strip single quotes around sheet name
            ws_name = ws_name[1:-1]
        region = wb[ws_name][reg]
    else:
        # passed a named range; find the cells in the workbook
        full_range = wb.get_named_range(range_name)
        if full_range is None:
            raise ValueError(
                'Range "{}" not found in workbook "{}".'.format(range_name, xlsx_file)
            )
        # convert to list (openpyxl 2.3 returns a list but 2.4+ returns a generator)
        destinations = list(full_range.destinations)
        if len(destinations) > 1:
            raise ValueError(
                'Range "{}" in workbook "{}" contains more than one region.'
                .format(range_name, xlsx_file)
            )
        ws, reg = destinations[0]
        # convert to worksheet object (openpyxl 2.3 returns a worksheet object
        # but 2.4+ returns the name of a worksheet)
        if isinstance(ws, str):
            ws = wb[ws]
        region = ws[reg]
    df = pd.DataFrame([cell.value for cell in row] for row in region)
    return df
