#!/usr/bin/env python
# coding: utf-8
"""This script contains the model formulation."""

# ############################ LICENSE INFORMATION ############################
# This file is part of the E3 RESOLVE Model.

# Copyright (C) 2019 Energy and Environmental Economics, Inc.
# For contact information, go to www.ethree.com

# The E3 RESOLVE Model is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.

# The E3 RESOLVE Model is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Affero General Public License for more details.

# You should have received a copy of the GNU Affero General Public License
# along with the E3 RESOLVE Model (in the file LICENSE.TXT). If not,
# see <http://www.gnu.org/licenses/>.
# #############################################################################

import logging
import pyomo.environ as pyo
from pyomo.core.base.PyomoModel import AbstractModel


class ResolveModel(AbstractModel):
    """RESOLVE model formulation and components.

    Includes intialization, validation, and fomrulation methods.

    Args:
        AbstractModel (pyo.AbstractModel)
    """

    def __init__(self, feature_toggles):
        """Initialize RESOLVE model with feature toggles.

        Args:
            feature_toggles (dict): Dictionary of feature settings defined in feature_toggles.csv
        """
        AbstractModel.__init__(self)

        # Assign feature toggle to ResolveModel
        for key, value in feature_toggles.items():
            setattr(self, key, value)

        # Create a custom logger
        self.logger = logging.getLogger('resolve_model.model_formulation')

        # Pyomo component calls
        self.construct_sets_and_params()
        self.construct_vars_and_expressions()
        self.construct_constraints()
        self.construct_objective()

    def construct_sets_and_params(self):
        """Construct sets and parameters."""
        ##################################################
        # TEMPORAL SETS AND PARAMETERS                   #
        ##################################################

        self.TIMEPOINTS = pyo.Set(domain=pyo.PositiveIntegers, ordered=True, doc="Unique study timepoints")
        self.period = pyo.Param(self.TIMEPOINTS, within=pyo.PositiveIntegers)
        self.month = pyo.Param(self.TIMEPOINTS, within=pyo.PositiveIntegers)
        self.day = pyo.Param(self.TIMEPOINTS, within=pyo.PositiveIntegers)
        self.hour_of_day = pyo.Param(self.TIMEPOINTS, within=pyo.NonNegativeIntegers)
        self.PERIODS = pyo.Set(domain=pyo.PositiveIntegers, initialize=_periods_init, ordered=True)
        self.VINTAGES = pyo.Set(domain=pyo.PositiveIntegers, initialize=_periods_init, ordered=True)
        self.first_period = pyo.Param(domain=pyo.PositiveIntegers, initialize=_first_period_init)
        self.PERIOD_VINTAGES = pyo.Set(
            dimen=2,
            within=self.PERIODS * self.VINTAGES,
            initialize=_period_vintage_init,
            ordered=True
        )
        self.DAYS = pyo.Set(domain=pyo.PositiveIntegers, initialize=_days_init, ordered=True, doc='Unique study days.')
        self.HOURS_OF_DAY = pyo.Set(domain=pyo.NonNegativeIntegers, initialize=_hours_of_day_init, ordered=True)
        self.first_timepoint_of_day = pyo.Param(self.PERIODS, self.DAYS, initialize=_first_timepoint_of_day_init)
        self.last_timepoint_of_day = pyo.Param(self.PERIODS, self.DAYS, initialize=_last_timepoint_of_day_init)
        self.previous_timepoint = pyo.Param(self.TIMEPOINTS, initialize=_previous_timepoint_init)
        self.next_timepoint = pyo.Param(self.TIMEPOINTS, initialize=_next_timepoint_init)
        self.MONTHS = pyo.Set(initialize=_months_init, ordered=True)

        self.discount_factor = pyo.Param(self.PERIODS, within=pyo.NonNegativeReals)
        self.years_in_period = pyo.Param(self.PERIODS, within=pyo.NonNegativeReals)
        self.day_weight = pyo.Param(self.DAYS, within=pyo.NonNegativeReals)
        self.hours_per_year = pyo.Param(initialize=8760.0)
        self.timepoints_per_day = pyo.Param(initialize=24)

        ##################################################
        # SPATIAL SETS AND PARAMETERS                    #
        ##################################################

        self.ZONES = pyo.Set(ordered=True)
        self.TRANSMISSION_LINES = pyo.Set(ordered=True)
        self.transmission_from = pyo.Param(self.TRANSMISSION_LINES, within=self.ZONES)
        self.transmission_to = pyo.Param(self.TRANSMISSION_LINES, within=self.ZONES)
        if self.allow_tx_reliability_contribution:
            self.TRANSMISSION_RELIABILITY_MAP = pyo.Set(dimen=3, ordered=True)
            self.transmission_direction = pyo.Param(self.TRANSMISSION_RELIABILITY_MAP, within=pyo.Integers)
            self.existing_tx_reliability_contribution = pyo.Param(self.TRANSMISSION_RELIABILITY_MAP, within=pyo.NonNegativeReals)
            self.new_tx_reliability_contribution = pyo.Param(
                self.TRANSMISSION_RELIABILITY_MAP,
                within=pyo.NonNegativeReals)
        self.SIMULTANEOUS_FLOW_GROUPS = pyo.Set(ordered=True)
        self.simultaneous_flow_limit_mw = pyo.Param(self.SIMULTANEOUS_FLOW_GROUPS, self.PERIODS, within=pyo.Reals)
        self.SIMULTANEOUS_FLOW_GROUP_LINES = pyo.Set(
            dimen=2,
            within=self.SIMULTANEOUS_FLOW_GROUPS * self.TRANSMISSION_LINES,
            ordered=True
        )
        self.direction = pyo.Param(self.SIMULTANEOUS_FLOW_GROUP_LINES, within=pyo.Integers)
        self.direction_abs = pyo.Param(
            self.SIMULTANEOUS_FLOW_GROUP_LINES,
            initialize=_direction_abs_init,
            validate=_direction_validate
        )

        ##################################################
        # RESOURCE SETS AND PARAMETERS                 #
        ##################################################

        self.RESOURCES = pyo.Set(ordered=True)
        self.technology = pyo.Param(self.RESOURCES, within=pyo.Any) # for results aggregation
        self.zone = pyo.Param(self.RESOURCES, within=self.ZONES)
        self.can_build_new = pyo.Param(self.RESOURCES, within=pyo.Boolean)
        self.capacity_limited = pyo.Param(self.RESOURCES, within=pyo.Boolean)
        self.can_retire = pyo.Param(self.RESOURCES, within=pyo.Boolean)
        self.thermal = pyo.Param(self.RESOURCES, within=pyo.Boolean)
        self.unit_commitment = pyo.Param(self.RESOURCES, within=pyo.Boolean)
        self.generate_at_max = pyo.Param(self.RESOURCES, within=pyo.Boolean)
        self.variable = pyo.Param(self.RESOURCES, within=pyo.Boolean)
        self.conventional_dr = pyo.Param(self.RESOURCES, within=pyo.Boolean)
        self.energy_efficiency = pyo.Param(self.RESOURCES, within=pyo.Boolean)
        self.storage = pyo.Param(self.RESOURCES, within=pyo.Boolean)
        self.increases_load = pyo.Param(self.RESOURCES, within=pyo.Boolean, validate=_increases_load_validate)
        self.load_only = pyo.Param(self.RESOURCES, within=pyo.Boolean, validate=_load_only_validate)
        self.enable_hourly_energy_tracking = pyo.Param(
            self.RESOURCES,
            within=pyo.Boolean,
            validate=_enable_hourly_energy_tracking_validate
        )

        # Resource types, distinguished by operational characteristics
        self.THERMAL_RESOURCES = pyo.Set(
            within=self.RESOURCES,
            initialize=self.RESOURCES,
            filter=(lambda model,
                           resource: model.thermal[resource]),
            validate=(
                lambda model,
                resource: model.unit_commitment[resource] +
                model.generate_at_max[resource] == 1
            ),
            ordered=True
        )
        self.UNIT_COMMITMENT_RESOURCES = pyo.Set(
            within=self.THERMAL_RESOURCES,
            initialize=self.THERMAL_RESOURCES,
            filter=(lambda model,
                    resource: model.unit_commitment[resource]),
            ordered=True
        )
        self.GENERATE_AT_MAX_RESOURCES = pyo.Set(
            within=self.RESOURCES,
            initialize=self.RESOURCES,
            filter=(lambda model,
                           resource: model.generate_at_max[resource]),
            ordered=True
        )

        self.CONVENTIONAL_DR_RESOURCES = pyo.Set(
            within=self.RESOURCES,
            initialize=self.RESOURCES,
            filter=(lambda model,
                           resource: model.conventional_dr[resource]),
            ordered=True
        )

        self.EE_PROGRAMS = pyo.Set(
            within=self.RESOURCES,
            initialize=self.RESOURCES,
            filter=(lambda model,
                           resource: model.energy_efficiency[resource]),
            ordered=True
        )

        self.RESOURCES_THAT_INCREASE_LOAD = pyo.Set(
            within=self.RESOURCES,
            initialize=self.RESOURCES,
            filter=(lambda model,
                           resource: model.increases_load[resource]),
            ordered=True
        )

        # Resources that can only increase demand (that only act as a load)
        self.LOAD_ONLY_RESOURCES = pyo.Set(
            within=self.RESOURCES_THAT_INCREASE_LOAD,
            initialize=self.RESOURCES_THAT_INCREASE_LOAD,
            filter=(lambda model,
                           resource: model.load_only[resource]),
            ordered=True
        )
        self.VARIABLE_RESOURCES = pyo.Set(
            within=self.RESOURCES,
            initialize=self.RESOURCES,
            filter=(lambda model,
                           resource: model.variable[resource]),
            ordered=True
        )
        self.STORAGE_RESOURCES = pyo.Set(
            within=self.RESOURCES,
            initialize=self.RESOURCES,
            filter=(lambda model,
                           resource: model.storage[resource]),
            ordered=True
        )

        self.RESOURCES_WITH_ENERGY_IN_STORAGE = pyo.Set(
            within=self.RESOURCES_THAT_INCREASE_LOAD,
            initialize=self.RESOURCES_THAT_INCREASE_LOAD,
            filter=(lambda model,
                           resource: model.enable_hourly_energy_tracking[resource]),
            ordered=True
        )

        # resources that Provide_Power_MW will be defined over
        self.RESOURCES_WITH_PROVIDE_POWER = (self.RESOURCES -
                                             self.LOAD_ONLY_RESOURCES -
                                             self.EE_PROGRAMS)

        self.RESOURCES_FOR_GENERIC_PMIN_AND_PMAX = (self.RESOURCES_WITH_PROVIDE_POWER
                                                    - self.VARIABLE_RESOURCES
                                                    - self.UNIT_COMMITMENT_RESOURCES)



        # New resources
        self.NEW_BUILD_RESOURCES = pyo.Set(
            within=self.RESOURCES,
            initialize=self.RESOURCES,
            filter=(lambda model,
                           resource: model.can_build_new[resource]),
            ordered=True
        )
        self.NEW_BUILD_STORAGE_RESOURCES = pyo.Set(
            within=self.NEW_BUILD_RESOURCES,
            initialize=self.NEW_BUILD_RESOURCES,
            filter=(lambda model,
                           resource: model.storage[resource]),
            ordered=True
        )

        # Capacity-limited resources
        self.CAPACITY_LIMITED_RESOURCES = pyo.Set(
            within=self.NEW_BUILD_RESOURCES,
            initialize=self.NEW_BUILD_RESOURCES,
            filter=(lambda model,
                           resource: model.capacity_limited[resource]),
            ordered=True
        )
        self.capacity_limit_mw = pyo.Param(self.CAPACITY_LIMITED_RESOURCES, self.PERIODS, within=pyo.NonNegativeReals)

        if self.custom_period_constraints:
            self.CUSTOM_CONSTRAINT_PERIOD_GROUP_MAP = pyo.Set(dimen=4, ordered=True)

            self.custom_constraint_period_multiplier = pyo.Param(
                self.CUSTOM_CONSTRAINT_PERIOD_GROUP_MAP,
                within=pyo.Reals
            )

            self.CUSTOM_CONSTRAINT_PERIOD_GROUPS = pyo.Set(
                dimen=1,
                initialize=custom_constraint_period_group_init,
                ordered=True
            )

            self.custom_constraint_period_limit = pyo.Param(
                self.CUSTOM_CONSTRAINT_PERIOD_GROUPS,
                self.PERIODS,
                default=float('inf'),
                doc='#flexibleparam'
            )

        if self.custom_timepoint_constraints:
            self.CUSTOM_CONSTRAINT_TIMEPOINT_GROUP_MAP = pyo.Set(dimen=4, ordered=True)

            self.custom_constraint_timepoint_multiplier = pyo.Param(
                self.CUSTOM_CONSTRAINT_TIMEPOINT_GROUP_MAP,
                within=pyo.Reals
            )

            self.CUSTOM_CONSTRAINT_TIMEPOINT_GROUPS = pyo.Set(
                dimen=1,
                initialize=custom_constraint_timepoint_group_init,
                ordered=True
            )

            self.custom_constraint_timepoint_limit = pyo.Param(
                self.CUSTOM_CONSTRAINT_TIMEPOINT_GROUPS,
                self.TIMEPOINTS,
                default=float('inf'),
                doc='#flexibleparam'
            )

        self.energy_budget_cap_factor = pyo.Param(
            self.RESOURCES_WITH_PROVIDE_POWER,
            self.DAYS,
            default='undefined',
            validate=_fraction_or_undefined_validate,
            doc='#flexibleparam',
            within=pyo.Any
        )

        self.min_gen_fraction = pyo.Param(
            self.RESOURCES_FOR_GENERIC_PMIN_AND_PMAX,
            self.TIMEPOINTS,
            default=0,
            within=pyo.PercentFraction,
            doc='#flexibleparam'
        )
        # max_provide_power_fraction limits power production to a fraction of the available capacity in each timepoint.
        # This limit is in addition to any limit imposed by available_capacity_fraction.
        # Limits imposed by available_capacity_fraction reflect a fraction of the resource being totally offline,
        # whereas limits imposed by max_provide_power_fraction further limit the power production from the online fraction.
        self.max_provide_power_fraction = pyo.Param(
            self.RESOURCES_FOR_GENERIC_PMIN_AND_PMAX, self.TIMEPOINTS,
            default=1,
            within=pyo.PercentFraction,
            doc='#flexibleparam'
        )
        self.can_produce_less_than_energy_budget = pyo.Param(
            self.RESOURCES_WITH_PROVIDE_POWER,
            default=False,
            within=pyo.Boolean,
            validate=_can_produce_less_than_energy_budget_validate,
            doc='#flexibleparam'
        )
        self.max_increase_load_fraction = pyo.Param(
            self.RESOURCES_THAT_INCREASE_LOAD, self.TIMEPOINTS,
            default=1,
            within=pyo.PercentFraction,
            doc='#flexibleparam'
        )

        self.min_increase_load_fraction = pyo.Param(
            self.RESOURCES_THAT_INCREASE_LOAD, self.TIMEPOINTS,
            default=0,
            within=pyo.PercentFraction,
            doc='#flexibleparam'
        )

        # this assumes at least 1-hour ramps will be constrained
        self.max_ramp_duration_to_constrain = pyo.Param(within=pyo.PositiveIntegers)

        self.RAMP_DURATIONS = pyo.Set(
            within=pyo.PositiveIntegers,
            initialize=_ramp_durations_init,
            ordered=True
        )

        self.ramp_up_intertimepoint_limit_fraction = pyo.Param(
            self.RESOURCES_WITH_PROVIDE_POWER,
            self.RAMP_DURATIONS,
            default='undefined',
            validate=_fraction_or_undefined_validate_no_resource_that_increase_load,
            within=pyo.Any
        )
        self.ramp_down_intertimepoint_limit_fraction = pyo.Param(
            self.RESOURCES_WITH_PROVIDE_POWER,
            self.RAMP_DURATIONS,
            default='undefined',
            validate=_fraction_or_undefined_validate_no_resource_that_increase_load,
            within=pyo.Any
        )

        # resources that have an energy budget cap factor defined
        # validate that if an energy budget is defined for one day, it is defined for all days
        self.RESOURCES_WITH_ENERGY_BUDGETS = pyo.Set(
            within=self.RESOURCES_WITH_PROVIDE_POWER,
            initialize=self.RESOURCES_WITH_PROVIDE_POWER,
            ordered=True,
            filter=(lambda model,
                           resource: any(
                model.energy_budget_cap_factor[resource, day] != "undefined" for day in model.DAYS)
                    ),
            validate=(lambda model,
                             resource: all(
                model.energy_budget_cap_factor[resource, day] != "undefined" for day in model.DAYS)
                      )
        )

        ##################################################
        # RESERVE RESOURCES                              #
        ##################################################

        self.RESERVE_RESOURCES = pyo.Set(within=(self.RESOURCES - self.VARIABLE_RESOURCES), ordered=True)
        self.can_provide_spin = pyo.Param(self.RESERVE_RESOURCES, within=pyo.Boolean)
        self.SPINNING_RESERVE_RESOURCES = pyo.Set(
            within=self.RESERVE_RESOURCES,
            initialize=self.RESERVE_RESOURCES,
            filter=(lambda model,
                           resource: model.can_provide_spin[resource]),
            ordered=True
        )
        self.can_provide_reg = pyo.Param(self.RESERVE_RESOURCES, within=pyo.Boolean)
        self.REGULATION_RESERVE_RESOURCES = pyo.Set(
            within=self.RESERVE_RESOURCES,
            initialize=self.RESERVE_RESOURCES,
            filter=(lambda model,
                           resource: model.can_provide_reg[resource]),
            ordered=True
        )
        self.can_provide_lf_reserves = pyo.Param(self.RESERVE_RESOURCES, within=pyo.Boolean)
        self.LOAD_FOLLOWING_RESERVE_RESOURCES = pyo.Set(
            within=self.RESERVE_RESOURCES,
            initialize=self.RESERVE_RESOURCES,
            filter=(lambda model,
                           resource: model.can_provide_lf_reserves[resource]),
            ordered=True
        )
        self.contributes_to_min_gen = pyo.Param(self.RESERVE_RESOURCES, within=pyo.Boolean)
        self.MINIMUM_GENERATION_RESOURCES = pyo.Set(
            within=self.RESERVE_RESOURCES,
            initialize=self.RESERVE_RESOURCES,
            filter=(lambda model,
                           resource: model.contributes_to_min_gen[resource]),
            ordered=True
        )
        self.contributes_to_freq_resp_total_req = pyo.Param(self.RESERVE_RESOURCES, within=pyo.Boolean)
        self.TOTAL_FREQ_RESP_RESOURCES = pyo.Set(
            within=self.RESERVE_RESOURCES,
            initialize=self.RESERVE_RESOURCES,
            filter=(lambda model,
                           resource: model.contributes_to_freq_resp_total_req[resource]),
            ordered=True
        )

        self.contributes_to_freq_resp_partial_req = pyo.Param(self.RESERVE_RESOURCES, within=pyo.Boolean)
        # All resources that contribute to the partial frequency response requirement must also contribute to the total
        self.PARTIAL_FREQ_RESP_RESOURCES = pyo.Set(
            within=self.RESERVE_RESOURCES,
            initialize=self.TOTAL_FREQ_RESP_RESOURCES,
            filter=(lambda model,
                           resource: model.contributes_to_freq_resp_partial_req[resource]),
            ordered=True
        )

        self.thermal_freq_response_fraction_of_commitment = pyo.Param(
            self.RESERVE_RESOURCES,
            within=pyo.PercentFraction
        )

        ##################################################
        # RELIABILITY RESOURCES                          #
        ##################################################

        if self.deliverability:
            # Pyomo DataPortal will read the whole `resource_deliverability_zones.tab` file as a set
            self.RESOURCE_DELIVERABILITY_MAP = pyo.Set(dimen=2)

            self.DELIVERABILITY_ZONES = pyo.Set(ordered=True, doc="deliverability zones",
                                                initialize=deliverability_zones_init
                                                )

        self.CAN_RETIRE_RESOURCES = pyo.Set(
            within=self.RESOURCES,
            initialize=self.RESOURCES,
            filter=(lambda model,
                           resource: model.can_retire[resource]),
            ordered=True,
            validate=can_retire_resources_validation_rule
        )
        self.CAN_RETIRE_RESOURCES_NEW = pyo.Set(
            within=self.CAN_RETIRE_RESOURCES & self.NEW_BUILD_RESOURCES,
            initialize=self.CAN_RETIRE_RESOURCES & self.NEW_BUILD_RESOURCES,
            ordered=True
        )

        self.RELIABILITY_ENERGY_LIMITED_RESOURCES = pyo.Set(
            initialize=self.STORAGE_RESOURCES | self.CONVENTIONAL_DR_RESOURCES,
            ordered=True)

        ##################################################
        # FUEL AND THERMAL PARAMETERS                    #
        ##################################################

        # Fuels set
        self.FUELS = pyo.Set(ordered=True)

        self.fuel = pyo.Param(self.THERMAL_RESOURCES, within=self.FUELS)
        self.fuel_burn_slope_mmbtu_per_mwh = pyo.Param(self.THERMAL_RESOURCES, within=pyo.NonNegativeReals)
        self.fuel_burn_intercept_mmbtu_per_hr = pyo.Param(self.THERMAL_RESOURCES, within=pyo.NonNegativeReals)

        # Differentiate between blend fuels and base fuels?
        self.is_blend_fuel = pyo.Param(self.FUELS, within=pyo.Boolean)

        # Define the base fuel and blend fuel set
        self.BLEND_FUELS = pyo.Set(initialize=self.FUELS, filter=(lambda model, fuel: model.is_blend_fuel[fuel]))
        self.BASE_FUELS = self.FUELS - self.BLEND_FUELS
        self.BLENDABLE_BASE_FUELS = pyo.Set(within = self.BASE_FUELS)

        self.BLEND_FUEL_RESOURCES = pyo.Set(
            within=self.THERMAL_RESOURCES,
            initialize=self.THERMAL_RESOURCES,
            filter=(lambda model, resource: model.fuel[resource] in model.BLENDABLE_BASE_FUELS),
            ordered=True
        )

        # define the blend fuel for each blendable base fuels
        if self.allow_fuel_blending:
            self.blend_fuel = pyo.Param(self.BLENDABLE_BASE_FUELS, within = self.BLEND_FUELS)



        ##################################################
        # UNIT_COMMITMENT PARAMETERS                     #
        ##################################################

        self.min_stable_level_fraction = pyo.Param(self.UNIT_COMMITMENT_RESOURCES, within=pyo.PercentFraction)
        self.ramp_rate_fraction = pyo.Param(self.UNIT_COMMITMENT_RESOURCES, within=pyo.NonNegativeReals)
        self.min_down_time_hours = pyo.Param(self.UNIT_COMMITMENT_RESOURCES,
                                             within=pyo.RangeSet(0,self.timepoints_per_day))
        self.min_up_time_hours = pyo.Param(self.UNIT_COMMITMENT_RESOURCES,
                                           within=pyo.RangeSet(0,self.timepoints_per_day))
        self.unit_size_mw = pyo.Param(self.UNIT_COMMITMENT_RESOURCES, within=pyo.NonNegativeReals)
        self.startup_cost_per_mw = pyo.Param(self.UNIT_COMMITMENT_RESOURCES, within=pyo.NonNegativeReals)
        self.shutdown_cost_per_mw = pyo.Param(self.UNIT_COMMITMENT_RESOURCES, within=pyo.NonNegativeReals)
        self.start_fuel_mmbtu_per_mw = pyo.Param(self.UNIT_COMMITMENT_RESOURCES, within=pyo.NonNegativeReals)
        self.avg_heat_rate_mmbtu_per_mwh = pyo.Param(self.THERMAL_RESOURCES, rule=define_avg_heat_rate)
        self.min_stable_level_heat_rate_mmbtu_per_mwh = pyo.Param(
            self.THERMAL_RESOURCES,
            rule=define_min_stable_level_heat_rate
        )
        self.commit_all_capacity = pyo.Param(
            self.UNIT_COMMITMENT_RESOURCES,
            self.TIMEPOINTS,
            within=pyo.Boolean,
            default=False,
            doc='#flexibleparam'
        )
        self.UNIT_COMMITMENT_RAMP_LIMITED_RESOURCES = pyo.Set(
            within=self.UNIT_COMMITMENT_RESOURCES,
            initialize=self.UNIT_COMMITMENT_RESOURCES,
            filter=(lambda model,
                           resource: model.ramp_rate_fraction[resource] < 1),
            ordered=True
        )

        ##################################################
        # DISPATCH PARAMETERS                            #
        ##################################################

        self.max_daily_capacity_factor = pyo.Param(
            self.RESOURCES_WITH_PROVIDE_POWER,
            self.PERIODS,
            self.DAYS,
            within=pyo.PercentFraction,
            default=1.0,
            doc='#flexibleparam'
        )

        self.min_daily_increase_load_cap_factor = pyo.Param(
            self.RESOURCES_THAT_INCREASE_LOAD,
            self.PERIODS,
            self.DAYS,
            within=pyo.PercentFraction,
            default=0.0,
            doc='#flexibleparam'
        )

        self.adjacency_h = pyo.Param(
            self.RESOURCES_THAT_INCREASE_LOAD & self.RESOURCES_WITH_PROVIDE_POWER,
            within=pyo.PositiveIntegers,
            default=9999,
            doc='Optional limit on how many hours before/after energy can be shifted. #flexibleparam',
            validate=_adjacency_h_validation_rule
        )

        # energy_brought_online_mwh and energy_taken_offline_mwh allow the user to add or subtract an exogenous amount of energy to the storage energy budget
        self.energy_taken_offline_mwh = pyo.Param(
            self.RESOURCES_WITH_ENERGY_IN_STORAGE,
            self.TIMEPOINTS,
            default=0,
            within=pyo.NonNegativeReals,
            doc='#flexibleparam'
        )

        self.energy_brought_online_mwh = pyo.Param(
            self.RESOURCES_WITH_ENERGY_IN_STORAGE,
            self.TIMEPOINTS,
            default=0,
            within=pyo.NonNegativeReals,
            doc='#flexibleparam'
        )

        ##################################################
        # CAPACITY & NEW BUILD PARAMETERS                #
        ##################################################

        self.planned_installed_capacity_mw = pyo.Param(
            self.RESOURCES,
            self.PERIODS,
            within=pyo.NonNegativeReals
        )
        self.min_operational_planned_capacity_mw = pyo.Param(
            self.RESOURCES,
            self.PERIODS,
            within=pyo.NonNegativeReals,
            validate=min_operational_planned_capacity_validate
        )
        self.planned_addition_mw = pyo.Param(
            self.RESOURCES,
            self.PERIODS,
            within=pyo.NonNegativeReals,
            initialize=_calculate_planned_capacity_addition_init
        )
        self.planned_subtraction_mw = pyo.Param(
            self.RESOURCES,
            self.PERIODS,
            within=pyo.NonNegativeReals,
            initialize=_calculate_planned_capacity_subtraction_init
        )
        self.cumulative_planned_capacity_additions_mw = pyo.Param(
            self.RESOURCES,
            self.PERIODS,
            within=pyo.NonNegativeReals,
            initialize=_cumulative_planned_capacity_additions_init
        )
        self.cumulative_planned_capacity_subtractions_mw = pyo.Param(
            self.RESOURCES,
            self.PERIODS,
            within=pyo.NonNegativeReals,
            initialize=_cumulative_planned_capacity_subtractions_init
        )
        # Fixed O&M costs for planned capacity. Can vary by period.
        self.planned_capacity_fixed_o_and_m_dollars_per_kw_yr = pyo.Param(
            self.RESOURCES,
            self.PERIODS,
            within=pyo.NonNegativeReals
        )
        self.min_cumulative_new_build_mw = pyo.Param(
            self.NEW_BUILD_RESOURCES,
            self.PERIODS,
            within=pyo.NonNegativeReals
        )
        # New resource annualized capital cost
        self.capital_cost_per_kw_yr = pyo.Param(self.NEW_BUILD_RESOURCES, self.VINTAGES, within=pyo.NonNegativeReals)
        # New resource annual Fixed O&M costs, which can vary by vintage
        self.new_capacity_fixed_o_and_m_dollars_per_kw_yr = pyo.Param(
            self.NEW_BUILD_RESOURCES,
            self.VINTAGES,
            within=pyo.NonNegativeReals
        )


        if self.deliverability:
            self.DELIVERABILITY_RESOURCES = pyo.Set(
                initialize=deliverability_resources_init,
                validate=deliverability_resources_validation_rule
            )

            self.RESOURCE_DELIVERABILITY_ZONES = pyo.Set(
                self.DELIVERABILITY_RESOURCES,
                within=self.DELIVERABILITY_ZONES,
                ordered=True,
                initialize=resource_deliverability_zones_init
            )

        ##################################################
        # RESOURCES WITH ENERGY IN STORAGE PARAMETERS    #
        ##################################################
        self.planned_storage_energy_capacity_mwh = pyo.Param(
            self.STORAGE_RESOURCES,
            self.PERIODS,
            within=pyo.NonNegativeReals
        )
        self.planned_storage_energy_capacity_fixed_o_and_m_dollars_per_kwh_yr = pyo.Param(
            self.STORAGE_RESOURCES,
            self.PERIODS,
            within=pyo.NonNegativeReals
        )
        self.planned_storage_charge_capacity_mw = pyo.Param(
            self.RESOURCES_WITH_ENERGY_IN_STORAGE,
            self.PERIODS,
            default=float("inf"),
            within=pyo.NonNegativeReals,
            doc='#flexibleparam'
        )
        self.min_energy_in_storage_in_timepoint_duration_h = pyo.Param(
            self.RESOURCES_WITH_ENERGY_IN_STORAGE,
            self.TIMEPOINTS,
            default=0.0,
            within=pyo.Reals,
            validate=min_energy_validation_rule,
            doc='#flexibleparam'
        )
        # In the future this can replace planned_storage_capacity_mwh and can be tagged to new storage resource with a fixed duration
        self.fixed_duration_h = pyo.Param(
            self.RESOURCES_WITH_ENERGY_IN_STORAGE - self.STORAGE_RESOURCES,
            default=9999,
            within=pyo.PositiveReals,
            doc='#flexibleparam'
        )
        # Define only for storage resources with new build allowed
        self.energy_storage_cost_per_kwh_yr = pyo.Param(
            self.NEW_BUILD_STORAGE_RESOURCES,
            self.VINTAGES,
            within=pyo.NonNegativeReals
        )
        self.energy_storage_charge_capacity_cost_per_kw_yr = pyo.Param(
            self.NEW_BUILD_STORAGE_RESOURCES,
            self.VINTAGES,
            within=pyo.NonNegativeReals,
            default=0,
            validate=_energy_storage_charge_capacity_cost_validate,
            doc='#flexibleparam'
        )

        self.new_energy_capacity_fixed_o_and_m_dollars_per_kwh_yr = pyo.Param(
            self.NEW_BUILD_STORAGE_RESOURCES,
            self.VINTAGES,
            within=pyo.NonNegativeReals
        )

        if self.deliverability:
            self.new_deliverability_cost_per_mw_yr = pyo.Param(
                self.DELIVERABILITY_ZONES,
                self.PERIODS,
                within=pyo.NonNegativeReals,
                default=0.,
                doc='#flexibleparam'
            )

        ##################################################
        # MAINTENANCE PARAMETERS                         #
        ##################################################
        # available_capacity_fraction is a way to knock out part of the generation/resource fleet in certain hours.
        # This means that in addition to the upper bound on provide power being decreased,
        # the energy budget of storage goes down, the aggregate PMin of a unit commitment resource goes down, etc.
        self.available_capacity_fraction = pyo.Param(
            self.RESOURCES,
            self.TIMEPOINTS,
            within=pyo.PercentFraction,
            default=1.0,
            doc='#flexibleparam'
        )

        ##################################################
        # VARIABLE COST                                  #
        ##################################################

        self.variable_cost_per_mwh = pyo.Param(
            self.RESOURCES,
            self.TIMEPOINTS,
            within=pyo.Reals,
            default=0.,
            doc='#flexibleparam'
        )

        self.variable_cost_increase_load_per_mwh = pyo.Param(
            self.RESOURCES_THAT_INCREASE_LOAD,
            self.TIMEPOINTS,
            within=pyo.Reals,
            default=0.,
            doc='#flexibleparam'
        )

        ##################################################
        # RENEWABLES PARAMETERS                          #
        ##################################################
        # shape is the hourly capacity factor of variable renewables
        self.shape = pyo.Param(self.VARIABLE_RESOURCES, self.DAYS, self.HOURS_OF_DAY, within=pyo.PercentFraction)

        # resource-specific limitations on the ability to curtail/dispatch variable renewables
        self.curtailable = pyo.Param(
            self.VARIABLE_RESOURCES,
            within=pyo.Boolean,
            default=1,
            doc='#flexibleparam'
        )

        # define the set of variable renewable resources that can be curtailed
        self.CURTAILABLE_VARIABLE_RESOURCES = pyo.Set(
            within=self.VARIABLE_RESOURCES,
            initialize=self.VARIABLE_RESOURCES,
            filter=(lambda model,
                           resource: model.curtailable[resource]),
            ordered=True
        )

        if self.itc_paired_storage:
            self.PAIRED_STORAGE_RESOURCES = pyo.Set(within=self.STORAGE_RESOURCES, ordered=True)

            # currently assume that storage is only paired with new build variable resources
            self.paired_supply_resource = pyo.Param(
                self.PAIRED_STORAGE_RESOURCES,
                within=(self.VARIABLE_RESOURCES & self.NEW_BUILD_RESOURCES)
            )

        ##################################################
        # STORAGE PARAMETERS                             #
        ##################################################

        self.charging_efficiency = pyo.Param(
            self.RESOURCES_WITH_ENERGY_IN_STORAGE,
            within=pyo.PercentFraction,
            default=1,
            doc='#flexibleparam'
        )
        self.discharging_efficiency = pyo.Param(
            self.RESOURCES_WITH_ENERGY_IN_STORAGE,
            within=pyo.PercentFraction,
            default=1,
            doc='#flexibleparam',
            validate=discharge_efficiency_validate
        )
        self.min_duration_h = pyo.Param(
            self.STORAGE_RESOURCES,
            within=pyo.NonNegativeReals,
            default=1,
            doc='#flexibleparam'
        )

        ##################################################
        # TRANSMISSION PARAMETERS                        #
        ##################################################

        self.min_flow_planned_mw = pyo.Param(self.TRANSMISSION_LINES, within=pyo.NonPositiveReals)
        self.max_flow_planned_mw = pyo.Param(self.TRANSMISSION_LINES, within=pyo.NonNegativeReals)
        self.ramp_constrained = pyo.Param(self.TRANSMISSION_LINES, within=pyo.Boolean)
        self.max_intertie_ramp_duration_to_constrain = pyo.Param(within=pyo.PositiveIntegers)
        self.new_build_tx_flag = pyo.Param(self.TRANSMISSION_LINES, within=pyo.Boolean)

        # Hurdle rate costs represent the cost to move power along each direction of a transmission line.
        # Must be non-negative.
        self.positive_direction_hurdle_rate_per_mw = pyo.Param(
            self.TRANSMISSION_LINES,
            self.PERIODS,
            within=pyo.NonNegativeReals
        )
        self.negative_direction_hurdle_rate_per_mw = pyo.Param(
            self.TRANSMISSION_LINES,
            self.PERIODS,
            within=pyo.NonNegativeReals
        )
        self.TRANSMISSION_LINES_NEW = pyo.Set(
            within=self.TRANSMISSION_LINES,
            initialize=self.TRANSMISSION_LINES,
            filter=(lambda model,
                           line: model.new_build_tx_flag[line]),
            ordered=True
        )
        self.RAMP_CONSTRAINED_TRANSMISSION_LINES = pyo.Set(
            within=self.TRANSMISSION_LINES,
            initialize=self.TRANSMISSION_LINES,
            filter=(lambda model,
                           line: model.ramp_constrained[line]),
            ordered=True
        )

        if self.transmission_ramp_limit:
            self.INTERTIE_FLOW_RAMP_DURATIONS = pyo.Set(
                within=pyo.PositiveIntegers,
                initialize=_flow_ramp_durations_init,
                ordered=True
            )
            self.flow_ramp_up_limit_fraction = pyo.Param(
                self.RAMP_CONSTRAINED_TRANSMISSION_LINES,
                self.INTERTIE_FLOW_RAMP_DURATIONS,
                within=pyo.PercentFraction
            )
            self.flow_ramp_down_limit_fraction = pyo.Param(
                self.RAMP_CONSTRAINED_TRANSMISSION_LINES,
                self.INTERTIE_FLOW_RAMP_DURATIONS,
                within=pyo.PercentFraction
            )


        # ## New tranmsmission line subset and params ## #
        if self.allow_tx_build:
            # Params defining minimum and maximum available build decisions for each period for new and existing paths
            self.max_tx_build_mw = pyo.Param(self.TRANSMISSION_LINES_NEW, self.PERIODS, within=pyo.NonNegativeReals)
            self.min_tx_build_mw = pyo.Param(self.TRANSMISSION_LINES_NEW, self.PERIODS, within=pyo.NonNegativeReals)
            # Cost params
            self.new_tx_fixed_cost_per_mw_yr = pyo.Param(
                self.TRANSMISSION_LINES_NEW,
                self.PERIODS,
                within=pyo.NonNegativeReals
            )
            self.new_tx_fixed_cost_annualized_per_yr = pyo.Param(
                self.TRANSMISSION_LINES_NEW,
                self.PERIODS,
                within=pyo.NonNegativeReals
            )


        ##################################################
        # CONVENTIONAL SHED DEMAND RESPONSE              #
        ##################################################

        self.conventional_dr_availability_hours_per_year = pyo.Param(
            self.CONVENTIONAL_DR_RESOURCES,
            self.PERIODS,
            within=pyo.NonNegativeReals
        )
        self.conventional_dr_reliability_hours_per_day = pyo.Param(
            self.CONVENTIONAL_DR_RESOURCES,
            self.PERIODS,
            within=pyo.NonNegativeReals
        )


        ##################################################
        # SYSTEM PARAMETERS                              #
        ##################################################

        self.input_load_mw = pyo.Param(self.ZONES, self.TIMEPOINTS, within=pyo.NonNegativeReals)
        self.unserved_energy_penalty_per_mw = pyo.Param(within=pyo.NonNegativeReals)
        self.allow_unserved_energy = pyo.Param(within=pyo.Boolean)
        self.overgen_penalty_per_mw = pyo.Param(within=pyo.NonNegativeReals)
        self.spin_violation_penalty_per_mw = pyo.Param(within=pyo.NonNegativeReals)
        self.upward_lf_reserve_violation_penalty_per_mw = pyo.Param(within=pyo.NonNegativeReals)
        self.downward_lf_reserve_violation_penalty_per_mw = pyo.Param(within=pyo.NonNegativeReals)
        self.upward_reg_violation_penalty_per_mw = pyo.Param(within=pyo.NonNegativeReals)
        self.downward_reg_violation_penalty_per_mw = pyo.Param(within=pyo.NonNegativeReals)
        self.reserve_timeframe_fraction_of_hour = pyo.Param(within=pyo.PercentFraction)

        ##################################################
        # RENEWABLE PORTFOLIO STANDARDS (RPS)            #
        ##################################################

        if self.include_RPS:
            self.RPS_POLICIES = pyo.Set(ordered = True)
            # retail sales
            self.retail_sales_mwh = pyo.Param(self.RPS_POLICIES, self.PERIODS, within=pyo.NonNegativeReals)
            # use to calculate rps target: RPS target as a percentage of retail sales
            self.rps_fraction_of_retail_sales = pyo.Param(self.RPS_POLICIES, self.PERIODS, within=pyo.NonNegativeReals)
            # non modeled rps energy that represents factors outside of the model
            self.rps_nonmodeled_mwh = pyo.Param(self.RPS_POLICIES, self.PERIODS, within=pyo.NonNegativeReals)
            # as a fraction of the RPS target, what is the limit on unbundled RECs?
            self.rps_unbundled_fraction_limit = pyo.Param(self.RPS_POLICIES, self.PERIODS, within=pyo.PercentFraction)

            # the fraction of generation from each resource that is eligible for each rps policy.
            # 0 = resource does not count towards the policy.  1 = all generation from the resource counts for the policy.
            # 0 < rps_eligible_fraction < 1 can be used to model partial ownership or other situations in which the output
            # of a resource cannot contribute fully to an RPS policy.
            self.rps_eligible_fraction = pyo.Param(self.RPS_POLICIES, self.RESOURCES, within=pyo.NonNegativeReals)

            # A set of rps eligible resources for each rps policy, defined by the rps_eligible parameter
            self.RPS_ELIGIBLE_RESOURCES = pyo.Set(
                self.RPS_POLICIES,
                within=self.RESOURCES,
                initialize= (lambda model, policy:
                            [r for r in model.RESOURCES if model.rps_eligible_fraction[policy, r]>0]),
                ordered=True
            )
            # EE and RPS_blend_fuels contribute to RPS programs in a special and separate way. The rest of the
            # RPS_Eligible_Reousrces is called RPS_REGULAR_RESOURCES
            self.RPS_REGULAR_RESOURCES = pyo.Set(
                self.RPS_POLICIES,
                within=self.RESOURCES,
                initialize= _rps_regular_resources_init,
                ordered=True
            )


            # A set covering all the resources that contribute to any of the RPS policies
            self.RPS_ELIGIBLE_RESOURCES_ALL = pyo.Set(
                within = self.RESOURCES,
                initialize = _rps_eligible_resources_all_init,
                ordered = True,
            )

            # As it stands, the load zones and RPS target has two interactions:
            # 1. Zone determine if a resource is counted as unbundled
            # 2. In zone storage losses mitigate RPS.
            self.is_zone_in_rps_policy = pyo.Param(self.RPS_POLICIES, self.ZONES, within=pyo.Boolean)

            # initialize RPS zones based on is_zone_in_rps_policy: The zone to RPS policy table
            self.RPS_ZONES = pyo.Set(
                self.RPS_POLICIES,
                within=self.ZONES,
                initialize=(lambda model, policy:
                            [z for z in model.ZONES if model.is_zone_in_rps_policy[policy, z]]),
                ordered=True
            )

            self.starting_rps_bank_mwh = pyo.Param(
                self.RPS_POLICIES,
                within=pyo.NonNegativeReals,
                validate=rps_starting_input_validate,
                doc='If optimized banking is allowed, define the starting amount of RPS bank.'
            )

            self.enforce_unbundled_fraction_limit = pyo.Param(self.RPS_POLICIES, self.PERIODS,
                                                              rule=define_unbundled_flag, within=pyo.Boolean)

        # cost of curtailment - the exogenously assumed cost at which different contract zones
        # would be willing to curtail their variable renewable generation
        self.curtailment_cost_per_mwh = pyo.Param(
            self.CURTAILABLE_VARIABLE_RESOURCES,
            self.PERIODS,
            within=pyo.NonNegativeReals,
            default=0,
            validate=curtailment_cost_validate,
            doc='#flexibleparam'
        )

        ##################################################
        # RESERVE PARAMETERS                             #
        ##################################################

        self.include_in_load_following = pyo.Param(self.ZONES, within=pyo.Boolean)
        self.LOAD_FOLLOWING_ZONES = pyo.Set(
            within=self.ZONES,
            initialize=self.ZONES,
            filter=(lambda model,
                           zone: model.include_in_load_following[zone]),
            ordered=True,
            doc='Zones for which the load is assumed to be covered under load following requirements'
        )

        self.LOAD_FOLLOWING_ZONE_CURTAILABLE_RESOURCES = pyo.Set(
            within=self.CURTAILABLE_VARIABLE_RESOURCES,
            initialize=self.CURTAILABLE_VARIABLE_RESOURCES,
            filter=(lambda model,
                           resource: model.zone[resource] in model.LOAD_FOLLOWING_ZONES)
        )

        def _spin_reserve_req_validate(model, value, timepoint):
            """Ensure that if spin reserve is entered by timepoint, the spin reserve by zone is zero"""
            if value != 0:
                return all(model.spin_reserve_fraction_of_load[zone] == 0 for zone in model.ZONES)
            else:
                return True

        self.spin_reserve_fraction_of_load = pyo.Param(self.ZONES, within=pyo.PercentFraction)
        self.spin_reserve_req_mw = pyo.Param(self.TIMEPOINTS, within=pyo.NonNegativeReals,
                                             validate=_spin_reserve_req_validate)
        self.upward_reg_req = pyo.Param(self.TIMEPOINTS, within=pyo.NonNegativeReals)
        self.downward_reg_req = pyo.Param(self.TIMEPOINTS, within=pyo.NonNegativeReals)
        self.upward_lf_reserve_req = pyo.Param(self.TIMEPOINTS, within=pyo.NonNegativeReals)
        self.downward_lf_reserve_req = pyo.Param(self.TIMEPOINTS, within=pyo.NonNegativeReals)
        self.min_gen_committed_mw = pyo.Param(self.TIMEPOINTS, within=pyo.NonNegativeReals)
        self.freq_resp_total_req_mw = pyo.Param(self.TIMEPOINTS, within=pyo.NonNegativeReals)
        self.freq_resp_partial_req_mw = pyo.Param(self.TIMEPOINTS, within=pyo.NonNegativeReals)

        self.resource_upward_lf_req = pyo.Param(
            self.VARIABLE_RESOURCES,
            self.DAYS,
            self.HOURS_OF_DAY,
            within=pyo.NonNegativeReals
        )

        self.resource_downward_lf_req = pyo.Param(
            self.VARIABLE_RESOURCES,
            self.DAYS,
            self.HOURS_OF_DAY,
            within=pyo.NonNegativeReals
        )

        # These parameters characterize subhourly reserve mileage,
        # the amount of energy dispatched within the hour per unit of reserve provision.
        self.reg_dispatch_fraction = pyo.Param(within=pyo.PercentFraction)
        self.lf_reserve_dispatch_fraction = pyo.Param(within=pyo.PercentFraction)

        # What fraction of curtailable variable renewable generation is available for load following reserves?
        # Currently variable renewables are not assumed to provide regulation
        self.var_rnw_available_for_lf_reserves = pyo.Param(within=pyo.PercentFraction)

        # What fraction of the downward load following reserve requirement can be met with renewables?
        self.max_var_rnw_lf_reserves = pyo.Param(within=pyo.PercentFraction)

        # ### Fuel prices ### #
        self.fuel_price_per_mmbtu = pyo.Param(self.FUELS, self.PERIODS, self.MONTHS, within=pyo.NonNegativeReals)

        ##################################################
        # TUNING PARAMETERS                              #
        ##################################################

        # Tuning params remove degeneracies or rounding errors in the solution
        # These should be small enough to not change decisions with non-zero cost
        # but large enough to not cause numerical problems
        self.ramp_relax = pyo.Param(initialize=10 ** -4)

        ##################################################
        # RELIABILITY CONSTRAINT PARAMETERS       #
        ##################################################

        self.RELIABILITY_CONSTRAINT_RESOURCE_MAP = pyo.Set(dimen=2, ordered=True)
        self.RELIABILITY_CONSTRAINTS = pyo.Set(
            dimen=1,
            initialize=lambda model: sorted(set([x[0] for x in model.RELIABILITY_CONSTRAINT_RESOURCE_MAP])),
            ordered=True)
        self.RELIABILITY_RESOURCES = pyo.Set(
            self.RELIABILITY_CONSTRAINTS,
            initialize=reliability_resources_init,
            ordered=True
        )
        # The fraction above peak load needed to ensure reliability. For example to model a 15% planning reserve margin, the reliability_margin is 0.15 (=15% above peak load).
        self.reliability_margin = pyo.Param(self.RELIABILITY_CONSTRAINTS, self.PERIODS, within=pyo.PercentFraction)
        self.reliability_peak_load_mw = pyo.Param(self.RELIABILITY_CONSTRAINTS, self.PERIODS, within=pyo.NonNegativeReals)
        self.reliability_capacity_adjustment_mw = pyo.Param(self.RELIABILITY_CONSTRAINTS, self.PERIODS, within=pyo.Reals)
        self.net_qualifying_capacity_fraction = pyo.Param(self.RELIABILITY_CONSTRAINT_RESOURCE_MAP,
                                                          self.PERIODS,
                                                          within=pyo.NonNegativeReals,
                                                          default=1,
                                                          validate=_nqc_validation_rule)
        self.reliability_event_length_hour = pyo.Param(self.RELIABILITY_CONSTRAINTS, self.PERIODS, within=pyo.PositiveReals)

        ##################################################
        # RENEWABLE ELCC SURFACE FACETS                  #
        ##################################################

        if self.allow_elcc_surfaces:
            # elcc surface mapping for reliability constraint, surface name and resource
            self.ELCC_SURFACE_MAP = pyo.Set(dimen=3, ordered=True)
            # axis index for each resource on the surface
            self.elcc_surface_axes = pyo.Param(
                self.ELCC_SURFACE_MAP,
                within=pyo.Integers
            )
            # multiplier for each entry in the elcc surface mapping
            self.elcc_surface_axes_map_multiplier = pyo.Param(
                self.ELCC_SURFACE_MAP,
                within=pyo.NonNegativeReals
            )
            # all elcc surface names
            self.ELCC_SURFACE_NAMES = pyo.Set(
                dimen=2,
                initialize=lambda model: sorted(set([x[:-1] for x in model.ELCC_SURFACE_MAP])),
                ordered=True
            )
            # all resources assigned to ELCC surfaces, used for final reliability accounting
            # ELCC resources would not count towards firm capacity in reliability calculation
            self.ELCC_RESOURCES = pyo.Set(
                self.RELIABILITY_CONSTRAINTS,
                initialize=elcc_resources_init,
                ordered=True
            )

            # count of axes for each surface name
            self.elcc_surface_axes_count = pyo.Param(
                self.ELCC_SURFACE_NAMES,
                initialize=elcc_surface_axes_count_rule,
                within=pyo.NonNegativeIntegers
            )
            # unique ELCC surface facets indexed by surface name, period, facet #, axis
            self.ELCC_SURFACE_FACETS_AXES_PERIOD = pyo.Set(
                dimen=5,
                ordered=True
            )

            # get surface_facet_axes, split period out
            self.ELCC_SURFACE_FACETS_AXES = pyo.Set(
                dimen=4,
                initialize=lambda model: sorted(set([x[:-1] for x in model.ELCC_SURFACE_FACETS_AXES_PERIOD])),
                ordered=True
            )

            self.ELCC_SURFACE_FACETS = pyo.Set(
                dimen=3,
                initialize=lambda model: sorted(set([x[:-2] for x in model.ELCC_SURFACE_FACETS_AXES_PERIOD])),
                ordered=True
            )

            # coefficient for each facet-axis, slopes for resources and intercepts
            self.elcc_surface_facets_axes_coefficient = pyo.Param(
                self.ELCC_SURFACE_FACETS_AXES_PERIOD,
                within = pyo.Reals
            )



        ##################################################
        # RELIABILITY CONSTRAINT                 #
        ##################################################

        # ##### Emissions target params ##### #

        if self.include_emissions:
            self.EMISSION_TYPES = pyo.Set(ordered=True)
            self.emission_type_units = pyo.Param(self.EMISSION_TYPES, within=pyo.Any)

            self.RESOURCE_EMISSIONS_GROUPS_MAP = pyo.Set(dimen=3, ordered=True)
            self.EMISSIONS_GROUPS = pyo.Set(
                dimen=1,
                initialize=emissions_groups_init,
                ordered=True
            )

            self.emissions_group_type = pyo.Param(
                self.EMISSIONS_GROUPS,
                initialize=emissions_group_type_init,
                within=pyo.Any
            )

            self.EMISSIONS_GROUP_CONSTRAINT_PERIODS = pyo.Set(dimen=2,ordered=True)

            if self.include_emissions_targets:
                self.emissions_target_units_per_year = pyo.Param(self.EMISSIONS_GROUP_CONSTRAINT_PERIODS, within=pyo.NonNegativeReals)

                # the emissions credit adjusts the target for factors outside of the modeled resources and transmission lines.
                # it can be positive (relaxing the target) or negative (making the target more stringent)
                self.emissions_credit_units_per_year = pyo.Param(self.EMISSIONS_GROUP_CONSTRAINT_PERIODS, within=pyo.Reals)

            self.resource_emission_factor = pyo.Param(self.RESOURCES_WITH_PROVIDE_POWER,
                                         self.EMISSION_TYPES,
                                         within=pyo.Reals,
                                         default=0)

            self.resource_emission_factor_blended_fuel = pyo.Param(self.RESOURCES_WITH_PROVIDE_POWER,
                                         self.EMISSION_TYPES,
                                         within=pyo.Reals,
                                         default=0)

            self.resource_emission_factor_denominator_mmbtu_or_mwh = pyo.Param(self.RESOURCES_WITH_PROVIDE_POWER,
                                                                  self.EMISSION_TYPES,
                                                                  default="mwh",
                                                                  validate=resource_emission_factor_denominator_validate,
                                                                  within=pyo.Any
                                                                  )

            self.resource_emission_factor_validate = pyo.Expression(self.EMISSIONS_GROUPS,
                                                       self.FUELS,
                                                       rule=resource_emission_factor_validate_rule
            )

            # emissions import factors representing the amount of emissions
            # assumed to be flowing along each direction of each transmission line.
            # Must be non-negative.
            self.emission_factor_positive_direction_per_mwh = pyo.Param(
                self.TRANSMISSION_LINES,
                self.EMISSIONS_GROUPS,
                self.PERIODS,
                within=pyo.NonNegativeReals,
                default=0
            )
            self.emission_factor_negative_direction_per_mwh = pyo.Param(
                self.TRANSMISSION_LINES,
                self.EMISSIONS_GROUPS,
                self.PERIODS,
                within=pyo.NonNegativeReals,
                default=0,
                validate=emission_factor_negative_direction_validate
            )


        ##################################################
        # MULTI-DAY ENERGY SHARING                       #
        ##################################################

        if self.interday_energy_sharing:
            self.interday_energy_sharing_interval_id = pyo.Param(self.DAYS, within=pyo.NonNegativeReals)

            self.ENERGY_SHARING_INTERVAL = pyo.Set(
                domain=pyo.NonNegativeIntegers,
                initialize=_energy_sharing_interval_init,
                ordered=True
            )

            self.daily_max_energy_budget_increase_hours = pyo.Param(
                self.RESOURCES_WITH_PROVIDE_POWER,
                self.DAYS,
                default=0,
                within=pyo.NonNegativeReals
            )

            self.daily_max_energy_budget_decrease_hours = pyo.Param(
                self.RESOURCES_WITH_PROVIDE_POWER,
                self.DAYS,
                default=0,
                within=pyo.NonNegativeReals
            )
            self.max_energy_to_share_hours = pyo.Param(
                self.RESOURCES_WITH_PROVIDE_POWER,
                self.ENERGY_SHARING_INTERVAL,
                default=0,
                within=pyo.NonNegativeReals
            )

            self.RESOURCES_WITH_MULTIDAY_SHARING = pyo.Set(
                within=self.RESOURCES_WITH_PROVIDE_POWER,
                initialize=self.RESOURCES_WITH_PROVIDE_POWER,
                ordered=True,
                filter=(lambda model,
                               resource: any(
                    model.max_energy_to_share_hours[resource, interval] > 0 for interval in
                    model.ENERGY_SHARING_INTERVAL
                )
                        )
            )

        ##################################################
        # ENERGY EFFICIENCY INVESTMENT                   #
        ##################################################

        if self.allow_ee_investment:
            self.ee_t_and_d_losses_fraction = pyo.Param(self.EE_PROGRAMS, within=pyo.PercentFraction)
            self.max_investment_in_period_aMW = pyo.Param(self.EE_PROGRAMS, self.VINTAGES, within=pyo.NonNegativeReals)
            self.ee_shapes_btm_mwh_per_amw = pyo.Param(
                self.EE_PROGRAMS,
                self.DAYS,
                self.HOURS_OF_DAY,
                within=pyo.NonNegativeReals
            )

        ##################################################
        # SEMI-STORAGE ZONES                             #
        ##################################################

        if self.allow_semi_storage_zones:
            self.SEMI_STORAGE_ZONES = pyo.Set(ordered=True)

            self.ssz_from_zone = pyo.Param(self.SEMI_STORAGE_ZONES)

            semi_storage_zones_params = [
                'ssz_positive_direction_hurdle_rate_per_mw',
                'ssz_negative_direction_hurdle_rate_per_mw'
            ]

            for param in semi_storage_zones_params:
                setattr(self, param, pyo.Param(self.SEMI_STORAGE_ZONES, self.PERIODS, within=pyo.NonNegativeReals))

            semi_storage_zones_reals_params = ['ssz_max_flow_mw', 'ssz_min_flow_mw']

            for param in semi_storage_zones_reals_params:
                setattr(self, param, pyo.Param(self.SEMI_STORAGE_ZONES, self.PERIODS, within=pyo.Reals))

        ##################################################
        #  DEDICATED RESOURCE TRANSMISSION USE           #
        ##################################################

        if self.resource_use_tx_capacity:
            self.RESOURCE_TX_IDS = pyo.Set(ordered=True)

            # the name of the dedicated import resources
            self.dedicated_import_resource = pyo.Param(self.RESOURCE_TX_IDS, within=pyo.Any)

            # the name the of tx line that the resources are used to transmit power
            self.tx_line_used = pyo.Param(self.RESOURCE_TX_IDS, within=pyo.Any)

            # which direction of transmission line is used by the resource (1 or -1)
            self.resource_tx_direction = pyo.Param(self.RESOURCE_TX_IDS, within=pyo.Integers)

    def construct_vars_and_expressions(self):
        """Construct variables and expressions."""

        ##################################################
        # BUILD AND RETIREMENT VARIABLES                 #
        ##################################################

        self.Build_Capacity_MW = pyo.Var(self.NEW_BUILD_RESOURCES, self.VINTAGES, within=pyo.NonNegativeReals)
        self.Build_Storage_Energy_Capacity_MWh = pyo.Var(
            self.NEW_BUILD_STORAGE_RESOURCES,
            self.VINTAGES,
            within=pyo.NonNegativeReals
        )
        self.Cumulative_New_Installed_Capacity_MW = pyo.Expression(
            self.NEW_BUILD_RESOURCES,
            self.PERIODS,
            rule=new_installed_capacity_tracking_rule
        )

        if self.allow_tx_reliability_contribution:
            self.Existing_Tx_Reliability_Contribution_MW = pyo.Expression(
                self.TRANSMISSION_RELIABILITY_MAP,
                rule=existing_tx_reliability_contribution_rule)

        # New storage charge capacity build
        self.New_Storage_Charge_Capacity_Build_MW = pyo.Var(
            self.NEW_BUILD_STORAGE_RESOURCES,
            self.PERIODS,
            within=pyo.NonNegativeReals)

        self.Cumulative_New_Storage_Charge_Capacity_Build_MW = pyo.Expression(
            self.NEW_BUILD_STORAGE_RESOURCES,
            self.PERIODS,
            rule=cumulative_new_storage_charge_capacity_build_tracking_rule
        )

        # Retirement variables for planned capacity
        self.Retire_Planned_Capacity_MW = pyo.Var(self.CAN_RETIRE_RESOURCES, self.PERIODS, within=pyo.NonNegativeReals)

        # Retirement variables for new capacity.
        # Fixed O&M costs can vary by vintage, so retirements must be tracked by vintage.
        # Represents all retirements of a given vintage made up to and including the period in question.
        self.Retire_New_Capacity_By_Vintage_Cumulative_MW = pyo.Var(
            self.CAN_RETIRE_RESOURCES_NEW,
            self.PERIOD_VINTAGES,
            within=pyo.NonNegativeReals
        )
        self.Retire_New_Capacity_In_Period_Cumulative_MW = pyo.Expression(
            self.CAN_RETIRE_RESOURCES_NEW,
            self.PERIODS,
            rule=retire_new_capacity_period_total_rule
        )
        self.Operational_New_Capacity_In_Period_By_Vintage_MW = pyo.Expression(
            self.NEW_BUILD_RESOURCES,
            self.PERIOD_VINTAGES,
            rule=operational_new_capacity_in_period_by_vintage_rule
        )
        self.Operational_New_Capacity_MW = pyo.Expression(
            self.NEW_BUILD_RESOURCES,
            self.PERIODS,
            rule=operational_new_capacity_in_period_rule
        )
        self.Retire_Planned_Capacity_Cumulative_MW = pyo.Expression(
            self.CAN_RETIRE_RESOURCES,
            self.PERIODS,
            rule=cumulative_retired_planned_capacity_rule
        )
        self.Operational_Planned_Capacity_MW = pyo.Expression(
            self.RESOURCES,
            self.PERIODS,
            rule=operational_planned_capacity_in_period_rule
        )
        self.Operational_Capacity_MW = pyo.Expression(
            self.RESOURCES,
            self.PERIODS,
            rule=operational_capacity_tracking_rule
        )
        self.Available_Capacity_In_Timepoint_MW = pyo.Expression(
            self.RESOURCES,
            self.TIMEPOINTS,
            rule=available_capacity_def
        )
        # specific to storage resources where charging capacity != operational capacity
        self.Operational_Charging_Capacity_MW = pyo.Expression(
            self.STORAGE_RESOURCES,
            self.PERIODS,
            rule=operational_charging_capacity_rule)
        self.Max_Increase_Load_In_Timepoint_MW = pyo.Expression(
            self.RESOURCES_THAT_INCREASE_LOAD,
            self.TIMEPOINTS,
            rule=max_increase_load_in_timepoint_rule
        )
        self.Min_Increase_Load_In_Timepoint_MW = pyo.Expression(
            self.RESOURCES_THAT_INCREASE_LOAD,
            self.TIMEPOINTS,
            rule=min_increase_load_in_timepoint_rule
        )
        self.Cumulative_New_Storage_Energy_Capacity_MWh = pyo.Expression(
            self.NEW_BUILD_STORAGE_RESOURCES,
            self.PERIODS,
            rule=new_storage_installed_energy_capacity_tracking_rule
        )
        self.Total_Storage_Energy_Capacity_MWh = pyo.Expression(
            self.STORAGE_RESOURCES,
            self.PERIODS,
            rule=total_storage_installed_energy_capacity_tracking_rule
        )

        if self.deliverability:
            self.Delivered_Capacity_MW = pyo.Var(
                self.DELIVERABILITY_RESOURCES,
                self.PERIODS,
                within=pyo.NonNegativeReals
            )
            self.Undelivered_Capacity_MW = pyo.Var(
                self.DELIVERABILITY_RESOURCES,
                self.PERIODS,
                within=pyo.NonNegativeReals
            )

        self.Reliability_Capacity_Delivered_MW = pyo.Expression(
                self.RESOURCES,
                self.PERIODS,
                rule=reliability_capacity_delivered_mw_rule)

        # Reliability Variables
        self.Reliability_Capacity_preNQC_MW = pyo.Var(self.RELIABILITY_CONSTRAINT_RESOURCE_MAP,
                                                      self.PERIODS,
                                                      within=pyo.NonNegativeReals)
        self.Reliability_Duration_MWh = pyo.Expression(self.RELIABILITY_ENERGY_LIMITED_RESOURCES,
                                                       self.PERIODS,
                                                       rule=reliability_duration_MWh_rule)

        # NQC defaults to 1 for everything
        # NQC_MW should apply to all resources
        self.Reliability_Capacity_NQC_MW = pyo.Expression(self.RELIABILITY_CONSTRAINT_RESOURCE_MAP,
                                                          self.PERIODS,
                                                          rule=reliability_capacity_nqc_derate_rule)

        if self.allow_elcc_surfaces:
            # ELCC variables
            self.ELCC_MW = pyo.Var(self.PERIODS, self.ELCC_SURFACE_NAMES, within=pyo.NonNegativeReals)
            # expression to export ELCC mw value for each facet, used in export_results.py
            self.Elcc_Surface_Facet_LHS = pyo.Expression(
                self.PERIODS,
                self.ELCC_SURFACE_FACETS,
                rule=elcc_surface_facet_LHS_rule
            )

        ##################################################
        # TRANSMISSION BUILD VARIABLES                   #
        ##################################################
        if self.allow_tx_build:
            # Build variables
            self.New_Tx_Period_Installed_Capacity_MW = pyo.Var(
                self.TRANSMISSION_LINES_NEW,
                self.PERIODS,
                within=pyo.NonNegativeReals
            )
            self.New_Tx_Period_Vintage_Installed_Capacity_MW = pyo.Expression(
                self.TRANSMISSION_LINES_NEW,
                self.PERIOD_VINTAGES,
                rule=new_tx_period_vintage_installed_capacity_mw
            )
            self.New_Tx_Total_Installed_Capacity_MW = pyo.Expression(
                self.TRANSMISSION_LINES_NEW,
                self.PERIODS,
                rule=new_tx_total_installed_capacity_mw
            )

            self.New_Tx_Max_Build = pyo.Constraint(
                self.TRANSMISSION_LINES_NEW,
                self.PERIODS,
                rule=new_tx_build_max_build
            )
            self.New_Tx_Min_Build = pyo.Constraint(
                self.TRANSMISSION_LINES_NEW,
                self.PERIODS,
                rule=new_tx_build_min_build
            )

            self.New_Tx_Period_Vintage_Build_Cost = pyo.Expression(
                self.TRANSMISSION_LINES_NEW,
                self.PERIOD_VINTAGES,
                rule=new_tx_build_cost
            )

            if self.allow_tx_reliability_contribution:
                self.New_Tx_Reliability_Contribution_MW = pyo.Expression(
                    self.TRANSMISSION_RELIABILITY_MAP,
                    rule=new_tx_reliability_contribution_rule)

        ##################################################
        # UNIT COMMITMENT AND DISPATCH VARIABLES         #
        ##################################################

        # Variable to track power provided by each resource
        # Will be restricted differently depending on resource tech type

        self.Provide_Power_MW = pyo.Var(
            self.RESOURCES_WITH_PROVIDE_POWER,
            self.TIMEPOINTS,
            within=pyo.NonNegativeReals
        )
        # Thermal units that are not must run (i.e. are unit commitment) have to be committed to provide power
        self.Commit_Units = pyo.Var(
            self.UNIT_COMMITMENT_RESOURCES,
            self.TIMEPOINTS,
            within=pyo.NonNegativeReals
        )  # may change to int
        self.Start_Units = pyo.Var(
            self.UNIT_COMMITMENT_RESOURCES,
            self.TIMEPOINTS,
            within=pyo.NonNegativeReals
        )  # may change to int
        self.Shut_Down_Units = pyo.Var(
            self.UNIT_COMMITMENT_RESOURCES,
            self.TIMEPOINTS,
            within=pyo.NonNegativeReals
        )  # may change to int
        self.PreStart_Units = pyo.Expression(
            self.UNIT_COMMITMENT_RESOURCES,
            self.TIMEPOINTS,
            rule=prestart_to_start_definition
        )
        self.PreShut_Down_Units = pyo.Expression(
            self.UNIT_COMMITMENT_RESOURCES,
            self.TIMEPOINTS,
            rule=preshut_down_to_shut_down_definition
        )
        self.Starting_Units = pyo.Expression(self.UNIT_COMMITMENT_RESOURCES, self.TIMEPOINTS, rule=define_starting_units)
        self.Shutting_Down_Units = pyo.Expression(
            self.UNIT_COMMITMENT_RESOURCES,
            self.TIMEPOINTS,
            rule=define_shutting_down_units
        )
        self.Fully_Operational_Units = pyo.Expression(
            self.UNIT_COMMITMENT_RAMP_LIMITED_RESOURCES,
            self.TIMEPOINTS,
            rule=fully_operational_units_def_rule
        )
        self.Start_Capacity_MW = pyo.Expression(
            self.UNIT_COMMITMENT_RESOURCES,
            self.TIMEPOINTS,
            rule=define_start_capacity
        )
        self.Commit_Capacity_MW = pyo.Expression(
            self.UNIT_COMMITMENT_RESOURCES,
            self.TIMEPOINTS,
            rule=define_commit_capacity
        )

        # Additional storage variables including charging of storage and variable to track available energy in storage
        # self.Charge_Storage_MW = pyo.Var(self.STORAGE_RESOURCES, self.TIMEPOINTS, within=pyo.NonNegativeReals)
        self.Energy_in_Storage_MWh = pyo.Var(
            self.RESOURCES_WITH_ENERGY_IN_STORAGE,
            self.TIMEPOINTS,
            within=pyo.Reals
        )
        self.Available_Energy_Capacity_MWh = pyo.Expression(
            self.RESOURCES_WITH_ENERGY_IN_STORAGE,
            self.TIMEPOINTS,
            rule=available_energy_capacity_mwh_rule
        )


        ##################################################
        # MULTI-DAY ENERGY SHARING VARIABLES       #
        ##################################################

        if self.interday_energy_sharing:
            self.Daily_Energy_Budget_Change_MWh = pyo.Var(
                self.RESOURCES_WITH_MULTIDAY_SHARING,
                self.PERIODS,
                self.DAYS,
                within=pyo.Reals,
                doc=
                'Daily energy budget increase can be positive and negative (positive: budget increase, negative: budget decrease)'
            )
            self.Positive_Energy_Budget_Moved_MWh = pyo.Var(
                self.RESOURCES_WITH_MULTIDAY_SHARING,
                self.PERIODS,
                self.DAYS,
                within=pyo.NonNegativeReals,
                doc='The absolute amount of energy budget being moved (=|Daily_Energy_Budget_Change_MWh|)'
            )

        ##################################################
        # ENERGY EFFICIENCY VARIABLES                    #
        ##################################################

        if self.allow_ee_investment:
            self.EE_Reduced_Load_FTM_MW = pyo.Expression(
                self.EE_PROGRAMS,
                self.TIMEPOINTS,
                rule=define_ee_program_load_in_timepoint
            )

        ##################################################
        # SEMI-STORAGE ZONES VARIABLES                   #
        ##################################################

        if self.allow_semi_storage_zones:
            # power transmit decision variables: positives means power flow from the "from zone" to the semi storage zone,
            # negative means power flow from the semi storage zone to the "from zone".
            # "from zones" are in the whole functional zone set. semi storage zone is a virtual concept act similarly as a
            # very large battery without efficiency losses but limit by charge and discharge capacity (MW). hurdle rates are
            # applied for transmitting between the from zones and ssz.
            # total power transmit need to be net 0 in the period
            self.SSZ_Transmit_Power_MW = pyo.Var(self.SEMI_STORAGE_ZONES, self.TIMEPOINTS, within=pyo.Reals)
            self.SSZ_Positive_Transmit_Power_MW = pyo.Var(
                self.SEMI_STORAGE_ZONES,
                self.TIMEPOINTS,
                within=pyo.NonNegativeReals
            )
            self.SSZ_Negative_Transmit_Power_MW = pyo.Var(
                self.SEMI_STORAGE_ZONES,
                self.TIMEPOINTS,
                within=pyo.NonNegativeReals
            )

        ##################################################
        # RESERVES VARIABLES                             #
        ##################################################

        ########### define upward reg

        def provide_upward_reg_net_mw_rule(model, resource, timepoint):
            """Define Upward Reg to include upward reg from charge, if storage."""
            if resource in model.STORAGE_RESOURCES:
                return model.Provide_Upward_Reg_MW[resource, timepoint] \
                       + model.Provide_Upward_Reg_From_Charge_MW[resource, timepoint]
            else:
                return model.Provide_Upward_Reg_MW[resource, timepoint]

        self.Provide_Upward_Reg_MW = pyo.Var(
            self.REGULATION_RESERVE_RESOURCES,
            self.TIMEPOINTS,
            within=pyo.NonNegativeReals
        )

        self.Provide_Upward_Reg_From_Charge_MW = pyo.Var(
            (self.REGULATION_RESERVE_RESOURCES & self.STORAGE_RESOURCES),
            self.TIMEPOINTS,
            within=pyo.NonNegativeReals
        )

        self.Provide_Upward_Reg_Net_MW = pyo.Expression(
            self.REGULATION_RESERVE_RESOURCES,
            self.TIMEPOINTS,
            rule=provide_upward_reg_net_mw_rule
        )

        #######

        ######## define downward reg

        def provide_downward_reg_net_mw_rule(model, resource, timepoint):
            """Define Downward Reg to include downward reg from charge, if storage."""
            if resource in model.STORAGE_RESOURCES:
                return model.Provide_Downward_Reg_MW[resource, timepoint] \
                       + model.Provide_Downward_Reg_From_Charge_MW[resource, timepoint]
            else:
                return model.Provide_Downward_Reg_MW[resource, timepoint]

        self.Provide_Downward_Reg_MW = pyo.Var(
            self.REGULATION_RESERVE_RESOURCES,
            self.TIMEPOINTS,
            within=pyo.NonNegativeReals
        )

        self.Provide_Downward_Reg_From_Charge_MW = pyo.Var(
            (self.REGULATION_RESERVE_RESOURCES & self.STORAGE_RESOURCES),
            self.TIMEPOINTS,
            within=pyo.NonNegativeReals
        )

        self.Provide_Downward_Reg_Net_MW = pyo.Expression(
            self.REGULATION_RESERVE_RESOURCES,
            self.TIMEPOINTS,
            rule=provide_downward_reg_net_mw_rule
        )

        #########

        ########### define LF upward reserve

        def provide_lf_upward_reserve_net_mw_rule(model, resource, timepoint):
            """Define LF upward reserves to include reserves from charge, if storage."""
            if resource in model.STORAGE_RESOURCES:
                return model.Provide_LF_Upward_Reserve_MW[resource, timepoint] \
                       + model.Provide_LF_Upward_Reserve_From_Charge_MW[resource, timepoint]
            else:
                return model.Provide_LF_Upward_Reserve_MW[resource, timepoint]

        self.Provide_LF_Upward_Reserve_MW = pyo.Var(
            self.LOAD_FOLLOWING_RESERVE_RESOURCES,
            self.TIMEPOINTS,
            within=pyo.NonNegativeReals
        )

        self.Provide_LF_Upward_Reserve_From_Charge_MW = pyo.Var(
            (self.LOAD_FOLLOWING_RESERVE_RESOURCES & self.STORAGE_RESOURCES),
            self.TIMEPOINTS,
            within=pyo.NonNegativeReals
        )

        self.Provide_LF_Upward_Reserve_Net_MW = pyo.Expression(
            self.LOAD_FOLLOWING_RESERVE_RESOURCES,
            self.TIMEPOINTS,
            rule=provide_lf_upward_reserve_net_mw_rule
        )

        #######

        ######## define LF downward reserve

        def provide_lf_downward_reserve_net_mw_rule(model, resource, timepoint):
            """Define LF downward reserves to include reserves from charge, if storage."""
            if resource in model.STORAGE_RESOURCES:
                return model.Provide_LF_Downward_Reserve_MW[resource, timepoint] \
                       + model.Provide_LF_Downward_Reserve_From_Charge_MW[resource, timepoint]
            else:
                return model.Provide_LF_Downward_Reserve_MW[resource, timepoint]

        self.Provide_LF_Downward_Reserve_MW = pyo.Var(
            self.LOAD_FOLLOWING_RESERVE_RESOURCES,
            self.TIMEPOINTS,
            within=pyo.NonNegativeReals
        )

        self.Provide_LF_Downward_Reserve_From_Charge_MW = pyo.Var(
            (self.LOAD_FOLLOWING_RESERVE_RESOURCES & self.STORAGE_RESOURCES),
            self.TIMEPOINTS,
            within=pyo.NonNegativeReals
        )

        self.Provide_LF_Downward_Reserve_Net_MW = pyo.Expression(
            self.LOAD_FOLLOWING_RESERVE_RESOURCES,
            self.TIMEPOINTS,
            rule=provide_lf_downward_reserve_net_mw_rule
        )

        #########

        ######## define spin

        def provide_spin_net_mw_rule(model, resource, timepoint):
            """Define Spin to include spin from charge, if storage."""
            if resource in model.STORAGE_RESOURCES:
                return model.Provide_Spin_MW[resource, timepoint] \
                       + model.Provide_Spin_From_Charge_MW[resource, timepoint]
            else:
                return model.Provide_Spin_MW[resource, timepoint]

        self.Provide_Spin_MW = pyo.Var(
            self.SPINNING_RESERVE_RESOURCES,
            self.TIMEPOINTS,
            within=pyo.NonNegativeReals
        )

        self.Provide_Spin_From_Charge_MW = pyo.Var(
            (self.SPINNING_RESERVE_RESOURCES & self.STORAGE_RESOURCES),
            self.TIMEPOINTS,
            within=pyo.NonNegativeReals
        )

        self.Provide_Spin_Net_MW = pyo.Expression(
            self.SPINNING_RESERVE_RESOURCES,
            self.TIMEPOINTS,
            rule=provide_spin_net_mw_rule
        )

        #########

        ######## define Freq. Resp.

        def provide_freq_resp_net_mw_rule(model, resource, timepoint):
            """Define freq. resp. to include freq. resp. from charge, if storage."""
            if resource in model.STORAGE_RESOURCES:
                return model.Provide_Frequency_Response_MW[resource, timepoint] \
                       + model.Provide_Frequency_Response_From_Charge_MW[resource, timepoint]
            else:
                return model.Provide_Frequency_Response_MW[resource, timepoint]

        self.Provide_Frequency_Response_MW = pyo.Var(
            self.TOTAL_FREQ_RESP_RESOURCES,
            self.TIMEPOINTS,
            within=pyo.NonNegativeReals
        )

        self.Provide_Frequency_Response_From_Charge_MW = pyo.Var(
            (self.TOTAL_FREQ_RESP_RESOURCES & self.STORAGE_RESOURCES),
            self.TIMEPOINTS,
            within=pyo.NonNegativeReals
        )

        self.Provide_Frequency_Response_Net_MW = pyo.Expression(
            self.TOTAL_FREQ_RESP_RESOURCES,
            self.TIMEPOINTS,
            rule=provide_freq_resp_net_mw_rule
        )
        #########


        # Violations
        self.Upward_Reg_Violation_MW = pyo.Var(self.TIMEPOINTS, within=pyo.NonNegativeReals)
        self.Downward_Reg_Violation_MW = pyo.Var(self.TIMEPOINTS, within=pyo.NonNegativeReals)
        self.Upward_LF_Reserve_Violation_MW = pyo.Var(self.TIMEPOINTS, within=pyo.NonNegativeReals)
        self.Downward_LF_Reserve_Violation_MW = pyo.Var(self.TIMEPOINTS, within=pyo.NonNegativeReals)
        self.Spin_Violation_MW = pyo.Var(self.TIMEPOINTS, within=pyo.NonNegativeReals)

        self.Scheduled_Curtailment_MW = pyo.Var(
            self.CURTAILABLE_VARIABLE_RESOURCES,
            self.TIMEPOINTS,
            within=pyo.NonNegativeReals,
            doc='Define an hourly (scheduled) curtailment variable for variable renewable resources that can be curtailed'
        )
        self.Variable_Resource_Provide_Downward_LF_MW = pyo.Var(
            self.TIMEPOINTS,
            within=pyo.NonNegativeReals,
            doc=
            'With endogenous load following, allow variable resources to provide upward load following (if pre-curtailed)'
        )
        self.Variable_Resource_Provide_Upward_LF_MW = pyo.Var(
            self.TIMEPOINTS,
            within=pyo.NonNegativeReals
        )
        self.Subhourly_Downward_LF_Energy_MWh = pyo.Expression(
            self.ZONES, self.TIMEPOINTS,
            rule=downward_lf_subhourly_energy_rule
        )
        self.Subhourly_Upward_LF_Energy_MWh = pyo.Expression(
            self.ZONES, self.TIMEPOINTS,
            rule=upward_lf_subhourly_energy_rule
        )
        self.Upward_Load_Following_Reserve_Req = pyo.Expression(
            self.TIMEPOINTS,
            rule=upward_combined_lf_reserve_req
        )
        self.Downward_Load_Following_Reserve_Req = pyo.Expression(
            self.TIMEPOINTS,
            rule=downward_combined_lf_reserve_req
        )

        ##################################################
        # TRANSMISSION VARIABLES                   #
        ##################################################

        self.Transmit_Power_Unspecified_MW = pyo.Var(self.TRANSMISSION_LINES, self.TIMEPOINTS, within=pyo.Reals)
        self.Transmit_Power_MW = pyo.Expression(self.TRANSMISSION_LINES, self.TIMEPOINTS, rule=define_transmit_power)

        # positive and negative direction version of Transmit_Power_Unspecified_MW
        # The constraints that define these variables assume that the variables are within pyo.NonNegativeReals,
        # which means that the negative direction variable will always be >= 0, even though
        # Transmit_Power_Unspecified_MW for the negative direction is negative (we're in effect taking ABS(Transmit_Power_MW)).

        # before adjusting for dedicated import resources that are located outside of the specified zone,
        # what is the positive and negative direction
        # this is used only in hurdle rates and emissions calculation: because the dedicated import resources are counted for emissions
        # and for hurdle rates directly, we assumed dedicated import resources have long term tx contract and we calculate
        # emissions costs directly
        self.Transmit_Power_Unspecified_Positive_Direction_MW = pyo.Var(
            self.TRANSMISSION_LINES,
            self.TIMEPOINTS,
            within=pyo.NonNegativeReals
        )
        self.Transmit_Power_Unspecified_Negative_Direction_MW = pyo.Var(
            self.TRANSMISSION_LINES,
            self.TIMEPOINTS,
            within=pyo.NonNegativeReals
        )
        self.Hurdle_Cost_Per_Timepoint_Positive_Direction = pyo.Expression(
            self.TRANSMISSION_LINES,
            self.TIMEPOINTS,
            rule=define_hurdle_cost_per_timepoint_positive_direction
        )
        self.Hurdle_Cost_Per_Timepoint_Negative_Direction = pyo.Expression(
            self.TRANSMISSION_LINES,
            self.TIMEPOINTS,
            rule=define_hurdle_cost_per_timepoint_negative_direction
        )

        if self.include_emissions:
            self.Unspecified_Imports_Units_Per_Timepoint_Positive_Direction = pyo.Expression(
                self.TRANSMISSION_LINES,
                self.EMISSIONS_GROUPS,
                self.TIMEPOINTS,
                rule=define_unspecified_import_emissions_per_timepoint_positive_direction
            )
            self.Unspecified_Imports_Units_Per_Timepoint_Negative_Direction = pyo.Expression(
                self.TRANSMISSION_LINES,
                self.EMISSIONS_GROUPS,
                self.TIMEPOINTS,
                rule=define_unspecified_import_emissions_per_timepoint_negative_direction
            )

            self.Unspecified_Imports_Units_Per_Year = pyo.Expression(
                self.EMISSIONS_GROUPS,
                self.PERIODS,
                rule=define_unspecified_import_emissions
            )

        ##################################################
        # DELIVERABILITY VARIABLES AND EXPRESSIONS       #
        ##################################################

        if self.deliverability:
            self.New_Deliverability_Incremental_MW = pyo.Var(self.DELIVERABILITY_ZONES, self.PERIODS,
                                                                  within=pyo.NonNegativeReals)
            self.New_Deliverability_Cumulative_MW = pyo.Expression(self.DELIVERABILITY_ZONES, self.PERIODS,
                                                                   rule=new_deliverability_capacity_tracking_rule)
            self.Deliverability_Cost_Annual_In_Period_Dollars = pyo.Expression(
                self.DELIVERABILITY_ZONES,
                self.PERIODS,
                rule=deliverability_cost_annual_in_period_dollars_rule)


        ##################################################
        # LOAD RESOURCE VARIABLES                        #
        ##################################################

        self.Increase_Load_MW = pyo.Var(
            self.RESOURCES_THAT_INCREASE_LOAD,
            self.TIMEPOINTS,
            within=pyo.NonNegativeReals
        )

        self.Zonal_Increase_Load_MW = pyo.Expression(
            self.ZONES,
            self.TIMEPOINTS,
            rule=sum_zonal_increase_load_rule
        )

        self.Zonal_Provide_Power_MW = pyo.Expression(
            self.ZONES,
            self.TIMEPOINTS,
            rule=sum_zonal_provide_power_rule
        )

        self.Zonal_Net_Provide_Power_MW = pyo.Expression(
            self.ZONES,
            self.TIMEPOINTS,
            rule=sum_zonal_net_provide_power_rule
        )

        ##################################################
        # SYSTEM VARIABLES                               #
        ##################################################

        self.Unserved_Energy_MW = pyo.Var(self.ZONES, self.TIMEPOINTS, within=pyo.NonNegativeReals)
        self.Overgeneration_MW = pyo.Var(self.ZONES, self.TIMEPOINTS, within=pyo.NonNegativeReals)

        ##################################################
        # FUEL VARIABLES                                 #
        ##################################################

        # Regular Fuel (Basis fuel as opposed to blend fuel)consumption
        self.Base_Fuel_Consumption_MMBtu = pyo.Var(
            self.THERMAL_RESOURCES,
            self.TIMEPOINTS,
            within=pyo.NonNegativeReals
        )
        # define fuel costs
        self.Base_Fuel_Cost = pyo.Expression(
            self.THERMAL_RESOURCES,
            self.TIMEPOINTS,
            rule=define_base_fuel_cost
        )

        if self.allow_fuel_blending:
            # Blend fuel consumption. Blend fuel are RPS eligible fuels that can be blended into fuel supply
            self.Blend_Fuel_Consumption_MMBtu = pyo.Var(
                self.BLEND_FUEL_RESOURCES,
                self.TIMEPOINTS,
                within=pyo.NonNegativeReals
            )
            self.Blend_Fuel_Cost = pyo.Expression(
                self.BLEND_FUEL_RESOURCES,
                self.TIMEPOINTS,
                rule=define_blend_fuel_cost
            )
            self.Blend_Fuel_Generation_MW = pyo.Expression(
                self.BLEND_FUEL_RESOURCES,
                self.TIMEPOINTS,
                rule = define_blend_fuel_generation_mw
            )

        self.Fuel_Consumption_MMBtu_Per_Year = pyo.Expression(
            self.FUELS, self.PERIODS,
            rule=define_yearly_fuel_consumption
        )

        self.Total_Fuel_Cost = pyo.Expression(
            self.THERMAL_RESOURCES,
            self.TIMEPOINTS,
            rule=define_total_fuel_cost
        )
        self.Start_Fuel_MMBtu = pyo.Expression(
            self.UNIT_COMMITMENT_RESOURCES,
            self.TIMEPOINTS,
            rule=define_start_fuel_in_timepoint
        )

        ##################################################
        # RPS VARIABLES                                  #
        ##################################################

        if self.include_RPS:
            self.Retail_Sales_After_EE_MWh = pyo.Expression(
                self.RPS_POLICIES,
                self.PERIODS,
                rule=define_sales_after_ee_mwh
            )
            self.RPS_Target_MWh = pyo.Expression(self.RPS_POLICIES, self.PERIODS, rule=define_rps_target)
            self.Bank_RPS_MWh = pyo.Var(self.RPS_POLICIES, self.PERIODS, within=pyo.NonNegativeReals)
            self.Previously_Banked_RPS_MWh = pyo.Expression(
                self.RPS_POLICIES, self.PERIODS,
                rule=define_previously_banked_rps)
            self.RPS_Storage_Losses_MWh_Per_Year = pyo.Expression(
                self.RPS_POLICIES, self.PERIODS,
                rule=define_yearly_rps_storage_losses
            )
            self.RPS_Eligible_Generation_MWh_Per_Year = pyo.Expression(
                self.RPS_POLICIES, self.PERIODS,
                rule=define_yearly_rps_eligible_generation
            )
            self.Scheduled_RPS_Curtailment_Per_Year = pyo.Expression(
                self.RPS_POLICIES, self.PERIODS,
                rule=define_scheduled_rps_curtailment_per_year
            )
            self.Subhourly_RPS_Curtailment_Per_Year = pyo.Expression(
                self.RPS_POLICIES, self.PERIODS,
                rule=define_subhourly_rps_curtailment_per_year
            )
            self.RPS_Blend_Fuel_Generation_MWh_Per_Year = pyo.Expression(
                self.RPS_POLICIES, self.PERIODS,
                rule=define_yearly_blend_fuel_rps_generation
            )
            self.Unbundled_RPS_Generation_MWh_Per_Year = pyo.Expression(
                self.RPS_POLICIES, self.PERIODS,
                rule = define_yearly_unbundled_rps_generation
            )

        ##################################################
        # EMISSIONS VARIABLES                                  #
        ##################################################

        if self.include_emissions:
            self.Resource_Emissions_Units_Per_Timepoint = pyo.Expression(
                self.RESOURCES_WITH_PROVIDE_POWER,
                self.EMISSION_TYPES,
                self.TIMEPOINTS,
                rule=define_resource_timepoint_emissions
            )

            self.Resource_Emissions_Units_Per_Year = pyo.Expression(self.EMISSIONS_GROUPS,
                                                                    self.PERIODS,
                                                                    rule=define_resource_emissions)

        ##################################################
        # COST EXPRESSIONS                               #
        ##################################################

        self.Variable_Cost_In_Timepoint = pyo.Expression(
            self.RESOURCES,
            self.TIMEPOINTS,
            rule=define_variable_cost_in_timepoint
        )
        self.Start_Cost_In_Timepoint = pyo.Expression(
            self.UNIT_COMMITMENT_RESOURCES,
            self.TIMEPOINTS,
            rule=define_start_cost_in_timepoint
        )
        self.Shutdown_Cost_In_Timepoint = pyo.Expression(
            self.UNIT_COMMITMENT_RESOURCES,
            self.TIMEPOINTS,
            rule=define_shutdown_cost_in_timepoint
        )

        self.Capital_Cost_Annual_In_Period_Dollars = pyo.Expression(
            self.NEW_BUILD_RESOURCES,
            self.PERIODS,
            rule=capital_cost_annual_in_period_dollars_rule
        )

        self.Fixed_OM_Annual_In_Period_Dollars = pyo.Expression(
            self.RESOURCES,
            self.PERIODS,
            rule=fixed_o_and_m_annual_in_period_dollars_rule
        )

        ##################################################
        # RELIABILITY CONSTRAINT EXPRESSIONS      #
        ##################################################

        if self.allow_tx_reliability_contribution:
            self.Transmission_Reliability_Contribution_MW = pyo.Expression(
                self.RELIABILITY_CONSTRAINTS,
                self.PERIODS,  # Add reliability constraint index once updated elsewhere
                rule=transmission_reliability_contribution_rule)

        self.Resource_Non_ELCC_Surface_Reliability_Contribution_MW = pyo.Expression(
            self.RELIABILITY_CONSTRAINTS,
            self.PERIODS,
            rule=resource_non_elcc_surface_reliability_contribution_rule
        )

        #########################################################
        # CUSTOM CONSTRAINTS: PERIOD and TIMEPOINT-based limits #
        #########################################################
        if self.custom_period_constraints:
            self.Custom_Constraint_Period_LHS = pyo.Expression(
                self.CUSTOM_CONSTRAINT_PERIOD_GROUPS,
                self.PERIODS,
                rule=custom_constraint_period_lhs
            )

        if self.custom_timepoint_constraints:
            self.Custom_Constraint_Timepoint_LHS = pyo.Expression(
                self.CUSTOM_CONSTRAINT_TIMEPOINT_GROUPS,
                self.TIMEPOINTS,
                rule=custom_constraint_timepoint_lhs
            )

    def construct_constraints(self):
        """Constructs all model constraints."""

        ##################################################
        # RELIABILITY AND ELCC CONSTRAINTS                       #
        ##################################################

        self.Reliability_Capacity_Constraint = pyo.Constraint(self.RELIABILITY_CONSTRAINTS, self.PERIODS, rule=reliability_capacity_period_rule)

        self.Reliability_preNQC_MW_Constraint = pyo.Constraint(
            self.RELIABILITY_CONSTRAINT_RESOURCE_MAP, self.PERIODS, rule=reliability_preNQC_power_rule)

        self.Reliability_preNQC_MWh_Constraint = pyo.Constraint(
            self.RELIABILITY_CONSTRAINT_RESOURCE_MAP, self.PERIODS, rule=reliability_preNQC_energy_rule)

        self.Reliability_preNQC_Conventional_DR_Annual_Constraint = pyo.Constraint(
            self.RELIABILITY_CONSTRAINT_RESOURCE_MAP,
            self.PERIODS,
            rule=reliability_preNQC_conventional_dr_annual_rule)

        if self.allow_elcc_surfaces:
            self.ELCC_Surface_Constraint = pyo.Constraint(self.PERIODS, self.ELCC_SURFACE_FACETS, rule=elcc_surface_rule)


        ##################################################
        # RETIREMENT CONSTRAINTS                         #
        ##################################################

        self.Minimum_Planned_Retirements_Constraint = pyo.Constraint(
            self.CAN_RETIRE_RESOURCES,
            self.PERIODS,
            rule=minimum_planned_retirements_rule
        )

        self.Maximum_Planned_Retirements_Constraint = pyo.Constraint(
            self.CAN_RETIRE_RESOURCES,
            self.PERIODS,
            rule=maximum_planned_retirements_rule
        )

        self.Retire_New_Increasing_Only_Constraint = pyo.Constraint(
            self.CAN_RETIRE_RESOURCES_NEW,
            self.PERIOD_VINTAGES,
            rule=retire_new_increasing_only_rule
        )

        self.Retire_New_Upper_Bound_Constraint = pyo.Constraint(
            self.CAN_RETIRE_RESOURCES_NEW,
            self.PERIOD_VINTAGES,
            rule=retire_new_upper_bound_rule
        )

        self.Minimum_Operational_Planned_Capacity_Constraint = pyo.Constraint(
            self.CAN_RETIRE_RESOURCES,
            self.PERIODS,
            rule=min_operational_planned_capacity_rule
        )

        self.Minimum_Cumulative_New_Build_Constraint = pyo.Constraint(
            self.NEW_BUILD_RESOURCES,
            self.PERIODS,
            rule=min_cumulative_new_build_rule
        )

        self.Capacity_Limit_Constraint = pyo.Constraint(
            self.CAPACITY_LIMITED_RESOURCES,
            self.PERIODS,
            rule=capacity_limit_rule
        )

        if self.custom_period_constraints:
            self.Custom_Constraint_Period = pyo.Constraint(
                self.CUSTOM_CONSTRAINT_PERIOD_GROUPS,
                self.PERIODS,
                rule=custom_constraint_period_rule
            )

        if self.custom_timepoint_constraints:
            self.Custom_Constraint_Timepoint = pyo.Constraint(
                self.CUSTOM_CONSTRAINT_TIMEPOINT_GROUPS,
                self.TIMEPOINTS,
                rule=custom_constraint_timepoint_rule
            )

        self.Min_Storage_Duration_Constraint = pyo.Constraint(
            self.NEW_BUILD_STORAGE_RESOURCES,
            self.PERIODS,
            rule=storage_minimum_duration_rule
        )

        ##################################################
        # FUEL CONSTRAINTS                               #
        ##################################################

        #The sum of regular fuel and blend fuel should equal to fuel requirement
        self.Fuel_Consumption_Constraint = pyo.Constraint(
            self.THERMAL_RESOURCES, self.TIMEPOINTS,
            rule= fuel_consumption_rule
        )

        ##################################################
        # RPS AND EMISSIONS CONSTRAINTS                        #
        ##################################################

        if self.include_RPS:
            self.Achieve_RPS_Constraint = pyo.Constraint(self.RPS_POLICIES, self.PERIODS, rule=achieve_rps_rule)

            self.RPS_Unbundled_Fraction_Limit_Constraint = pyo.Constraint(
                self.RPS_POLICIES, self.PERIODS,
                rule=rps_unbundled_fraction_limit_rule
            )

        if self.include_emissions & self.include_emissions_targets:
            self.Emissions_Target_Constraint = pyo.Constraint(self.EMISSIONS_GROUP_CONSTRAINT_PERIODS,
                                                              rule=emissions_target_rule)

        ##################################################
        # DISPATCH CONSTRAINTS                           #
        ##################################################

        self.Renewable_Power_Constraint = pyo.Constraint(
            self.VARIABLE_RESOURCES,
            self.TIMEPOINTS,
            rule=variable_power_rule
        )

        self.Maximum_Thermal_Commitment_Constraint = pyo.Constraint(
            self.UNIT_COMMITMENT_RESOURCES,
            self.TIMEPOINTS,
            rule=unit_commitment_commitment_rule
        )

        self.Unit_Commitment_Max_Gen_Up_Reserve_Constraint = pyo.Constraint(
            self.UNIT_COMMITMENT_RESOURCES,
            self.TIMEPOINTS,
            rule=unit_commitment_max_gen_and_up_reserves_rule
        )

        self.Unit_Commitment_Upward_Reserve_Ramp_Constraint = pyo.Constraint(
            self.UNIT_COMMITMENT_RESOURCES,
            self.TIMEPOINTS,
            rule=unit_commitment_upward_reserve_ramp_rule
        )

        self.Unit_Commitment_Downward_Reserve_Ramp_Constraint = pyo.Constraint(
            self.UNIT_COMMITMENT_RESOURCES,
            self.TIMEPOINTS,
            rule=unit_commitment_downward_reserve_ramp_rule
        )
        # Constraint applied over intersection of unit commitment resources and resources that can provide frequency response
        self.Unit_Commitment_Freq_Resp_Limit_Constraint = pyo.Constraint(
            self.UNIT_COMMITMENT_RESOURCES & self.TOTAL_FREQ_RESP_RESOURCES,
            self.TIMEPOINTS,
            rule=unit_commitment_freq_resp_limit
        )

        self.Thermal_Min_Gen_Down_Reserve_Constraint = pyo.Constraint(
            self.UNIT_COMMITMENT_RESOURCES,
            self.TIMEPOINTS,
            rule=unit_commitment_min_gen_and_down_reserves_rule
        )

        self.Unit_Commitment_Start_Constraint = pyo.Constraint(
            self.UNIT_COMMITMENT_RESOURCES,
            self.TIMEPOINTS,
            rule=unit_commitment_start_rule
        )

        self.Unit_Commitment_Shut_Down_Constraint = pyo.Constraint(
            self.UNIT_COMMITMENT_RESOURCES,
            self.TIMEPOINTS,
            rule=unit_commitment_shut_down_rule
        )

        self.Commitment_Tracking_Constraint = pyo.Constraint(
            self.UNIT_COMMITMENT_RESOURCES,
            self.TIMEPOINTS,
            rule=unit_commitment_resource_commitment_tracking_rule
        )

        self.Unit_Commitment_Resource_Ramp_Up_Constraint = pyo.Constraint(
            self.UNIT_COMMITMENT_RAMP_LIMITED_RESOURCES,
            self.TIMEPOINTS,
            rule=unit_commitment_resource_ramp_up_rule
        )

        self.Unit_Commitment_Resource_Ramp_Down_Constraint = pyo.Constraint(
            self.UNIT_COMMITMENT_RAMP_LIMITED_RESOURCES,
            self.TIMEPOINTS,
            rule=unit_commitment_resource_ramp_down_rule
        )

        self.Generate_At_Max_Constraint = pyo.Constraint(
            self.GENERATE_AT_MAX_RESOURCES,
            self.TIMEPOINTS,
            rule=generate_at_max
        )

        self.Main_Energy_Budget_Constraint = pyo.Constraint(
            (self.RESOURCES_WITH_ENERGY_BUDGETS - self.STORAGE_RESOURCES),
            self.PERIODS,
            self.DAYS,
            rule=energy_budget_rule
        )

        self.Max_Daily_Provide_Power_Constraint = pyo.Constraint(
            self.RESOURCES_WITH_PROVIDE_POWER,
            self.PERIODS,
            self.DAYS,
            rule=max_daily_provide_power_rule
        )

        self.Non_Unit_Commitment_Max_Gen_Up_Reserve_Constraint = pyo.Constraint(
            self.RESOURCES_FOR_GENERIC_PMIN_AND_PMAX,
            self.TIMEPOINTS,
            rule=non_unit_commitment_max_gen_and_up_reserves_rule
        )

        self.Non_Unit_Commitment_Min_Gen_Down_Reserve_Constraint = pyo.Constraint(
            self.RESOURCES_FOR_GENERIC_PMIN_AND_PMAX,
            self.TIMEPOINTS,
            rule=non_unit_commitment_min_gen_and_down_reserves_rule
        )

        self.Min_Increase_Load_And_Up_Reserves_Constraint = pyo.Constraint(
            self.RESOURCES_THAT_INCREASE_LOAD,
            self.TIMEPOINTS,
            rule=min_increase_load_and_up_reserves_rule
        )

        self.Max_Increase_Load_And_Down_Reserves_Constraint = pyo.Constraint(
            self.RESOURCES_THAT_INCREASE_LOAD,
            self.TIMEPOINTS,
            rule=max_increase_load_and_down_reserves_rule
        )

        # Rule for resources that shift energy around with a time adjacency limit
        self.Energy_Shift_Increase_Load_Adjacency_Constraint = pyo.Constraint(
            self.RESOURCES_THAT_INCREASE_LOAD & self.RESOURCES_WITH_PROVIDE_POWER,
            self.TIMEPOINTS,
            rule=energy_shift_increase_load_adjacency_rule
        )
        self.Energy_Shift_Provide_Power_Adjacency_Constraint = pyo.Constraint(
            self.RESOURCES_THAT_INCREASE_LOAD & self.RESOURCES_WITH_PROVIDE_POWER,
            self.TIMEPOINTS,
            rule=energy_shift_provide_power_adjacency_rule
        )

        self.Intertimepoint_Ramp_Constraint_LB = pyo.Constraint(
            self.RESOURCES_WITH_PROVIDE_POWER,
            self.TIMEPOINTS,
            self.RAMP_DURATIONS,
            rule=intertimepoint_ramp_rule_lb
        )

        self.Intertimepoint_Ramp_Constraint_UB = pyo.Constraint(
            self.RESOURCES_WITH_PROVIDE_POWER,
            self.TIMEPOINTS,
            self.RAMP_DURATIONS,
            rule=intertimepoint_ramp_rule_ub
        )

        self.Storage_Energy_Upper_Constraint = pyo.Constraint(
            self.RESOURCES_WITH_ENERGY_IN_STORAGE,
            self.TIMEPOINTS,
            rule=storage_energy_upper_rule
        )

        self.Storage_Energy_Lower_Constraint = pyo.Constraint(
            self.RESOURCES_WITH_ENERGY_IN_STORAGE,
            self.TIMEPOINTS,
            rule=storage_energy_lower_rule
        )

        self.Storage_Energy_Tracking_Constraint = pyo.Constraint(
            self.RESOURCES_WITH_ENERGY_IN_STORAGE,
            self.TIMEPOINTS,
            rule=storage_energy_tracking_rule
        )




        if self.itc_paired_storage:
            self.Paired_Storage_Charging_Constraint = pyo.Constraint(
                self.PAIRED_STORAGE_RESOURCES,
                self.TIMEPOINTS,
                rule=paired_battery_charging
            )

        ##################################################
        # POWER BALANCE CONSTRAINT                       #
        ##################################################

        self.Zonal_Power_Balance_Constraint = pyo.Constraint(self.ZONES, self.TIMEPOINTS, rule=zonal_power_balance_rule)

        ##################################################
        # TRANSMISSION CONSTRAINTS                       #
        ##################################################

        self.Transmission_Min_Flow_Constraint = pyo.Constraint(
            self.TRANSMISSION_LINES,
            self.TIMEPOINTS,
            rule=transmission_min_flow_rule
        )

        self.Transmission_Max_Flow_Constraint = pyo.Constraint(
            self.TRANSMISSION_LINES,
            self.TIMEPOINTS,
            rule=transmission_max_flow_rule
        )

        if self.transmission_ramp_limit:
            self.Transmission_Ramp_Down_Constraint = pyo.Constraint(
                self.RAMP_CONSTRAINED_TRANSMISSION_LINES,
                self.TIMEPOINTS,
                self.INTERTIE_FLOW_RAMP_DURATIONS,
                rule=transmission_ramp_down_rule
            )

            self.Transmission_Ramp_Up_Constraint = pyo.Constraint(
                self.RAMP_CONSTRAINED_TRANSMISSION_LINES,
                self.TIMEPOINTS,
                self.INTERTIE_FLOW_RAMP_DURATIONS,
                rule=transmission_ramp_up_rule
            )

        self.Simultaneous_Flows_Limit_Constraint = pyo.Constraint(
            self.SIMULTANEOUS_FLOW_GROUPS,
            self.TIMEPOINTS,
            rule=simultaneous_transmission_flows_rule
        )

        ##################################################
        # OPERATIONAL RESERVES CONSTRAINTS               #
        ##################################################

        self.Spinning_Reserve_Req_MW = pyo.Expression(self.TIMEPOINTS, rule=spinning_reserve_req)

        self.Meet_Spin_Requirement_Constraint = pyo.Constraint(self.TIMEPOINTS, rule=spinning_reserve_rule)

        self.Meet_Upward_Reg_Requirement_Constraint = pyo.Constraint(self.TIMEPOINTS, rule=upward_regulation_rule)

        self.Meet_Downward_Reg_Requirement_Constraint = pyo.Constraint(self.TIMEPOINTS, rule=downward_regulation_rule)

        self.Meet_Upward_LF_Requirement_Constraint = pyo.Constraint(
            self.TIMEPOINTS,
            rule=meet_upward_load_following_reserve_rule
        )

        self.Meet_Downward_LF_Requirement_Constraint = pyo.Constraint(
            self.TIMEPOINTS,
            rule=downward_load_following_reserve_rule
        )

        self.Variable_Resource_Available_Upward_LF_Constraint = pyo.Constraint(
            self.TIMEPOINTS,
            rule=upward_load_following_reserve_rule
        )

        self.Var_Renw_Down_LF_Reserve_Availability_Constraint = pyo.Constraint(
            self.TIMEPOINTS,
            rule=variable_rnw_down_reserve_availability_rule
        )

        self.Max_Upward_LF_From_Variable_Resources_Constraint = pyo.Constraint(
            self.TIMEPOINTS,
            rule=max_upward_lf_from_variable_rule
        )

        self.Var_Renw_Down_LF_Reserve_Limit_Constraint = pyo.Constraint(
            self.TIMEPOINTS,
            rule=variable_rnw_downward_lf_reserve_limit_rule
        )

        self.Total_Frequency_Response_Headroom_Constraint = pyo.Constraint(
            self.TIMEPOINTS,
            rule=total_frequency_response_rule
        )

        self.Partial_Frequency_Response_Headroom_Constraint = pyo.Constraint(
            self.TIMEPOINTS,
            rule=partial_frequency_response_rule
        )

        self.Min_Local_Gen_Constraint = pyo.Constraint(self.TIMEPOINTS, rule=minimum_local_committed_generation_rule)

        self.Storage_Discharge_Reserves_Energy_Constraint = pyo.Constraint(
            self.RESOURCES_WITH_ENERGY_IN_STORAGE,
            self.TIMEPOINTS,
            rule=storage_discharge_reserves_energy_rule
        )

        self.Storage_Charge_Reserves_Energy_Constraint = pyo.Constraint(
            self.RESOURCES_WITH_ENERGY_IN_STORAGE,
            self.TIMEPOINTS,
            rule=storage_charge_reserves_energy_rule
        )

        ##################################################
        # DELIVERABILITY CONSTRAINTS                     #
        ##################################################

        if self.deliverability:
            self.Deliverability_Definition_Constraint = pyo.Constraint(
                self.DELIVERABILITY_RESOURCES,
                self.PERIODS,
                rule=deliverability_def_rule
            )

            self.Delivered_Capacity_Period_Build_Constraint = pyo.Constraint(
                self.DELIVERABILITY_RESOURCES,
                self.PERIODS,
                rule=delivered_capacity_period_build_rule
            )

            self.Undelivered_Capacity_Period_Build_Constraint = pyo.Constraint(
                self.DELIVERABILITY_RESOURCES,
                self.PERIODS,
                rule=undelivered_capacity_period_build_rule
            )

        ##################################################
        # TRANSMISSION FLOW CONSTRAINTS                  #
        ##################################################

        self.Transmit_Power_Positive_Direction_Constraint = pyo.Constraint(
            self.TRANSMISSION_LINES,
            self.TIMEPOINTS,
            rule=positive_direction_transmit_power_rule
        )

        self.Transmit_Power_Negative_Direction_Constraint = pyo.Constraint(
            self.TRANSMISSION_LINES,
            self.TIMEPOINTS,
            rule=negative_direction_transmit_power_rule
        )

        ##################################################
        # DEMAND SIDE RESOURCES                          #
        ##################################################

        self.Min_Daily_Increase_Load_Constraint = pyo.Constraint(
            self.RESOURCES_THAT_INCREASE_LOAD,
            self.PERIODS,
            self.DAYS,
            rule=min_daily_increase_load_rule
        )

        self.Conventional_DR_Annual_Availability_Constraint = pyo.Constraint(
            self.CONVENTIONAL_DR_RESOURCES,
            self.PERIODS,
            rule=conventional_dr_max_annual_availability_rule
        )


        if self.allow_ee_investment:
            self.Maximum_EE_Investment_In_Period_Constraint = pyo.Constraint(
                self.EE_PROGRAMS,
                self.PERIODS,
                rule=maximum_ee_investment_in_period_rule
            )

        ##################################################
        # MULTI-DAY ENERGY SHARING CONSTRAINTS     #
        ##################################################

        if self.interday_energy_sharing:
            self.Net_Zero_Energy_Sharing_Constraint = pyo.Constraint(
                self.RESOURCES_WITH_MULTIDAY_SHARING,
                self.PERIODS,
                self.ENERGY_SHARING_INTERVAL,
                rule=net_zero_energy_sharing_rule
            )

            self.Min_Daily_Changes_Constraint = pyo.Constraint(
                self.RESOURCES_WITH_MULTIDAY_SHARING,
                self.PERIODS,
                self.DAYS,
                rule=min_daily_changes_rule
            )

            self.Max_Daily_Changes_Constraint = pyo.Constraint(
                self.RESOURCES_WITH_MULTIDAY_SHARING,
                self.PERIODS,
                self.DAYS,
                rule=max_daily_changes_rule
            )

            self.Define_Absolute_Energy_Moved_Constraint = pyo.Constraint(
                self.RESOURCES_WITH_MULTIDAY_SHARING,
                self.PERIODS,
                self.DAYS,
                rule=define_absolute_energy_moved_rule
            )

            self.Max_Absolute_Energy_Moved_Constraints = pyo.Constraint(
                self.RESOURCES_WITH_MULTIDAY_SHARING,
                self.PERIODS,
                self.ENERGY_SHARING_INTERVAL,
                rule=max_absolute_energy_moved_rule
            )

        ##################################################
        # SEMI-STORAGE ZONES CONSTRAINTS                 #
        ##################################################

        if self.allow_semi_storage_zones:
            self.SSZ_Energy_Net_0_Constraint = pyo.Constraint(
                self.SEMI_STORAGE_ZONES,
                self.PERIODS,
                rule=ssz_energy_net_0_rule
            )
            self.SSZ_Transmit_Power_Positive_Definition = pyo.Constraint(
                self.SEMI_STORAGE_ZONES,
                self.TIMEPOINTS,
                rule=ssz_transmit_power_positive_definition_rule
            )
            self.SSZ_Transmit_Power_Negative_Definition = pyo.Constraint(
                self.SEMI_STORAGE_ZONES,
                self.TIMEPOINTS,
                rule=ssz_transmit_power_negative_definition_rule
            )
            self.SSZ_Tx_Limit_Max_Constraint = pyo.Constraint(
                self.SEMI_STORAGE_ZONES,
                self.TIMEPOINTS,
                rule=ssz_tx_limit_max_rule
            )
            self.SSZ_Tx_Limit_Min_Constraint = pyo.Constraint(
                self.SEMI_STORAGE_ZONES,
                self.TIMEPOINTS,
                rule=ssz_tx_limit_min_rule
            )

    def construct_objective(self):
        """Constructs the objective function, which minimizes the NPV total cost of investment and operations."""
        self.Total_Cost = pyo.Objective(rule=total_cost_rule, sense=pyo.minimize)


##################################################
# COMPONENT DEFINITION HELPER METHODS            #
##################################################


# Investment periods
def _periods_init(model):
    """

    Args:
      model:

    Returns:

    """
    periods = list()
    for tmp_p in model.TIMEPOINTS:
        periods.append(model.period[tmp_p])

    return sorted(set(periods))


# The set PERIODS and the set VINTAGES are the same
# but VINTAGES denotes when a resource is built (as opposed to operated)


# find the first period
def _first_period_init(model):
    """

    Args:
      model:

    Returns:

    """
    return min(model.PERIODS)


def find_prev_period(model, period):
    """Returns the previous period, or none if the current period is the first period

    Args:
      model:
      period:

    Returns:

    """
    previous_periods = []
    for q in model.PERIODS:
        if q < period:
            previous_periods.append(q)

    if previous_periods != []:
        prev_period = max(previous_periods)
    else:
        prev_period = None

    return prev_period


def _period_vintage_init(model):
    """The set of tuples with the first element being the period and the second
    element being all vintages smaller than the period (i.e. capacity installed before the current period)

    Args:
      model:

    Returns:

    """
    period_vintages = list()
    for p in model.PERIODS:
        for v in model.VINTAGES:
            if v <= p:
                period_vintages.append((p, v))
    return period_vintages


def _days_init(model):
    """Days -- unique IDs within each period

    Args:
      model:

    Returns:

    """
    days = list()
    for tmp_d in model.TIMEPOINTS:
        days.append(model.day[tmp_d])

    return sorted(set(days))


def _hours_of_day_init(model):
    """Hours of day -- unique IDs within each day

    Args:
      model:

    Returns:

    """
    hours_of_day = list()
    for tmp in model.TIMEPOINTS:
        hours_of_day.append(model.hour_of_day[tmp])

    return sorted(set(hours_of_day))


def _first_timepoint_of_day_init(model, period, day):
    """Defines the first timepoint on each day.  Assumes timepoints are ordered.

    Args:
      model:
      period:
      day:
      period:

    Returns:

    """
    timepoints_on_day = list()
    for timepoint in model.TIMEPOINTS:
        if model.period[timepoint] == period and model.day[timepoint] == day:
            timepoints_on_day.append(timepoint)
    return min(timepoints_on_day)


def _last_timepoint_of_day_init(model, period, day):
    """Defines the last timepoint on each day.  Assumes timepoints are ordered.

    Args:
      model:
      day:
      period:

    Returns:

    """
    timepoints_on_day = list()
    for timepoint in model.TIMEPOINTS:
        if model.period[timepoint] == period and model.day[timepoint] == day:
            timepoints_on_day.append(timepoint)
    return max(timepoints_on_day)


def _previous_timepoint_init(model):
    """Define a "previous timepoint" for periodic boundary constraints
    The previous timepoint for the first hour of the day is the last hour of that day

    Args:
      model:

    Returns:

    """
    previous_tp = dict()
    for timepoint in model.TIMEPOINTS:
        if timepoint == model.first_timepoint_of_day[model.period[timepoint], model.day[timepoint]]:
            previous_tp[timepoint] = model.last_timepoint_of_day[model.period[timepoint], model.day[timepoint]]
        else:
            previous_tp[timepoint] = timepoint - 1
    return previous_tp


def _next_timepoint_init(model):
    """Define a "next timepoint" for periodic boundary constraints
    The next timepoint for the last hour of the day is the first hour of that day

    Args:
      model:

    Returns:

    """
    next_tp = dict()
    for timepoint in model.TIMEPOINTS:
        if timepoint == model.last_timepoint_of_day[model.period[timepoint], model.day[timepoint]]:
            next_tp[timepoint] = model.first_timepoint_of_day[model.period[timepoint], model.day[timepoint]]
        else:
            next_tp[timepoint] = timepoint + 1
    return next_tp


def _months_init(model):
    """Initialize months.

    Args:
      model:

    Returns:

    """
    months = []
    for tmp_m in model.TIMEPOINTS:
        months.append(model.month[tmp_m])

    return sorted(set(months))


# Make sure direction is -1 or 1
def _direction_abs_init(model, sfg, line):
    """

    Args:
      model:
      sfg:
      line:

    Returns:

    """
    return abs(model.direction[(sfg, line)])


def _direction_validate(model, a, sfg, line):
    """

    Args:
      model:
      a:
      sfg:
      line:

    Returns:

    """
    return a == 1


# create a list of unique groups from the larger map set
def custom_constraint_period_group_init(model):
    return sorted(set(
        val[0] for val in model.CUSTOM_CONSTRAINT_PERIOD_GROUP_MAP
    ))


def custom_constraint_timepoint_group_init(model):
    return sorted(set(
        val[0] for val in model.CUSTOM_CONSTRAINT_TIMEPOINT_GROUP_MAP
    ))

# Validation rule for resources that can be retired
# Rule is necessary because subsequent code can handle retirements for some resources but not others
# Storage resources can't be retired because doing so would require retiring energy and power capacity
# and the current code only retires power capacity.
def can_retire_resources_validation_rule(model, resource):
    """

    Args:
      model:
      resource:

    Returns:

    """
    if resource in model.THERMAL_RESOURCES - model.STORAGE_RESOURCES:
        return True
    else:
        return False


def resource_deliverability_zones_init(model, deliverability_resource):
    return sorted(set(
        val[1]
        for val in model.RESOURCE_DELIVERABILITY_MAP
        if val[0] == deliverability_resource
    ))


def deliverability_resources_init(model):
    return sorted(set(
        val[0]
        for val in model.RESOURCE_DELIVERABILITY_MAP
    ))

def deliverability_resources_validation_rule(model, resource):
    """ Deliverability functionality currently includes only resources that
        A) Do not have planned capacity (planned_installed_capacity_mw = 0 for all periods)
        B) Cannot retire (not in CAN_RETIRE_RESOURCES)
        C) Can be newly built (in NEW_BUILD_RESOURCES)
    """
    if (resource in model.NEW_BUILD_RESOURCES
        and resource not in model.CAN_RETIRE_RESOURCES
        and not any(model.planned_installed_capacity_mw[resource, period] != 0 for period in model.PERIODS)):
        return True
    else:
        return False


def deliverability_zones_init(model):
    return sorted(set(
        val[1]
        for val in model.RESOURCE_DELIVERABILITY_MAP
    ))


def reliability_resources_init(model, constraint):
    resource_list = []
    for val in model.RELIABILITY_CONSTRAINT_RESOURCE_MAP:
        if val[0]==constraint:
            resource_list.append(val[1])
    return resource_list

##############################################################################
# ####################### CALCULATED PARAMETERS ############################ #
##############################################################################

# ##### Heat rate params ##### #


def define_avg_heat_rate(model, resource):
    """Calculate the theoretical full load heat rate for each thermal resource.
    Due to the pseudo-MILP nature of RESOLVE, this is also the in-model heat rate during most operation conditions.
    Currently used for blend fuel and (in results reporting) production cost inputs

    Args:
      model:
      resource:

    Returns:

    """

    if resource in model.UNIT_COMMITMENT_RESOURCES:
        return ((model.fuel_burn_intercept_mmbtu_per_hr[resource]
                 / model.unit_size_mw[resource])
                + model.fuel_burn_slope_mmbtu_per_mwh[resource])
    else:
        return model.fuel_burn_slope_mmbtu_per_mwh[resource]


def define_min_stable_level_heat_rate(model, resource):
    """Calculate the min stable level heat rate for each thermal tech.
    Currently used for blend fuel and (in results reporting) production cost inputs

    Args:
      model:
      resource:

    Returns:

    """

    if resource in model.UNIT_COMMITMENT_RESOURCES:
        min_stable_level = model.min_stable_level_fraction[resource]
        if min_stable_level == 0:
            min_stable_level_heat_rate = 0
        else:
            min_capacity_mw = min_stable_level * model.unit_size_mw[resource]

            min_stable_level_heat_rate = \
                (model.fuel_burn_intercept_mmbtu_per_hr[resource]
                 + min_capacity_mw
                 * model.fuel_burn_slope_mmbtu_per_mwh[resource]) \
                / min_capacity_mw
    else:
        min_stable_level_heat_rate = model.fuel_burn_slope_mmbtu_per_mwh[resource]

    return min_stable_level_heat_rate


# ##### Resource params ##### #

def min_operational_planned_capacity_validate(model, min_operational_planned_capacity, resource, period):
    """If a resource can't be retired then make sure that there is a zero value for min_operational_planned_capacity_mw
    because any value would not be meaningful
    If a resource can be retired, make sure that the minimum amount of planned capacity
    that the model must keep operational in a period is less than or equal to the planned capacity

    Args:
      model:
      min_operational_planned_capacity:
      resource:
      period:

    Returns:

    """
    if resource not in model.CAN_RETIRE_RESOURCES:
        return min_operational_planned_capacity == 0
    else:
        return min_operational_planned_capacity <= model.planned_installed_capacity_mw[resource, period]


# Calculate planned additions for resources
def _calculate_planned_capacity_addition_init(model, resource, period):
    """

    Args:
      model:
      resource:
      period:

    Returns:

    """
    previous_period = find_prev_period(model, period)

    if period == model.first_period:
        capacity_addition = model.planned_installed_capacity_mw[resource, period]
    elif model.planned_installed_capacity_mw[resource, period] > \
            model.planned_installed_capacity_mw[resource, previous_period]:
        capacity_addition = model.planned_installed_capacity_mw[resource, period] - \
                            model.planned_installed_capacity_mw[resource, previous_period]
    else:
        capacity_addition = 0

    return capacity_addition


# Calculate planned subtractions (planned retirements) for resources.
# If the model can retire the resource, at least this amount of capacity will be retired before or in the period
def _calculate_planned_capacity_subtraction_init(model, resource, period):
    """

    Args:
      model:
      resource:
      period:

    Returns:

    """
    previous_period = find_prev_period(model, period)

    if period == model.first_period:
        capacity_subtraction = 0
    elif model.planned_installed_capacity_mw[resource, period] < \
            model.planned_installed_capacity_mw[resource, previous_period]:
        capacity_subtraction = model.planned_installed_capacity_mw[resource, previous_period] - \
                               model.planned_installed_capacity_mw[resource, period]
    else:
        capacity_subtraction = 0

    return capacity_subtraction


# Parameters tracking cumulative planned capacity additions for every period
def _cumulative_planned_capacity_additions_init(model, resource, period):
    """

    Args:
      model:
      resource:
      period:

    Returns:

    """
    cumulative_planned_capacity_additions = float()
    for p in model.PERIODS:
        if p <= period:
            cumulative_planned_capacity_additions += model.planned_addition_mw[resource, p]

    return cumulative_planned_capacity_additions


# Parameters tracking cumulative planned capacity subtractions (forced retirements) for every period
def _cumulative_planned_capacity_subtractions_init(model, resource, period):
    """

    Args:
      model:
      resource:
      period:

    Returns:

    """
    cumulative_planned_capacity_subtractions = float()
    for p in model.PERIODS:
        if p <= period:
            cumulative_planned_capacity_subtractions += model.planned_subtraction_mw[resource, p]

    return cumulative_planned_capacity_subtractions


def _ramp_durations_init(model):
    """
    Initialize ramp durations
    Args:
      model:

    Returns:

    """
    durations = list()
    for duration in range(1, model.max_ramp_duration_to_constrain + 1):
        durations.append(duration)
    return durations


def _flow_ramp_durations_init(model):
    """

    Args:
      model:

    Returns:

    """
    durations = list()
    for duration in range(1, model.max_intertie_ramp_duration_to_constrain + 1):
        durations.append(duration)
    return durations


def new_tx_period_vintage_installed_capacity_mw(model, line, period, vintage):
    """

    Args:
      model:
      line:
      period:
      vintage:

    Returns:

    """
    return model.New_Tx_Period_Installed_Capacity_MW[line, vintage]


def new_tx_total_installed_capacity_mw(model, line, period):
    """

    Args:
      model:
      line:
      period:

    Returns:

    """
    return sum(
        model.New_Tx_Period_Installed_Capacity_MW[line,
                                                  vintage] for vintage in model.VINTAGES if vintage <= period
    )


# Build constraints
def new_tx_build_max_build(model, line, period):
    """

    Args:
      model:
      line:
      period:

    Returns:

    """
    return (model.New_Tx_Total_Installed_Capacity_MW[line, period] <= model.max_tx_build_mw[line, period])


def new_tx_build_min_build(model, line, period):
    """

    Args:
      model:
      line:
      period:

    Returns:

    """
    return (model.New_Tx_Total_Installed_Capacity_MW[line, period] >= model.min_tx_build_mw[line, period])


# Build cost expression
def new_tx_build_cost(model, l, p, v):
    """

    Args:
      model:
      l:
      p:
      v:

    Returns:

    """
    return model.new_tx_fixed_cost_per_mw_yr[l, v] * model.New_Tx_Period_Vintage_Installed_Capacity_MW[l, p, v]


# ##### System params ##### #

# ### RPS ### #

def _rps_regular_resources_init(model, policy):
    ''' EE and blend_fuel can contribute to RPS, but their accounting is different and done somewhere else
    Create a set of RPS resources that are regular, i.e. excluding EE and blend_fuel

    :param model:
    :param policy:
    :return: the set of regular resources of each policy
    '''
    rps_regular_resources = []
    for r in model.RPS_ELIGIBLE_RESOURCES[policy]:
        if r not in (model.EE_PROGRAMS|model.BLEND_FUEL_RESOURCES):
            rps_regular_resources.append(r)

    return rps_regular_resources

def _rps_eligible_resources_all_init(model):
    """ Initialize a list of all the eligibible resources covered by at least one RPS.
    These resources are not allowed to have a curtailment cost. Only used in the curtailment cost validation

    :param model:
    :return rps_eligible_resources_all:
    """

    rps_eligible_resources_all = []
    for resource in model.RESOURCES:
        for rps_policy in model.RPS_POLICIES:
            if resource in model.RPS_ELIGIBLE_RESOURCES[rps_policy]:
                rps_eligible_resources_all.append(resource)
                break

    return rps_eligible_resources_all


def rps_starting_input_validate(model, starting_bank, policy):
    """The starting_rps_bank_mwh param is only used with optimized banking
    so check that this param is zero unless optimize_rps_banking is turned on

    Args:
      model:
      starting_bank:
      policy:

    Returns:

    """
    if not model.optimize_rps_banking:
        return starting_bank == 0
    else:
        return starting_bank >= 0


def curtailment_cost_validate(model, curtailment_cost, resource, period):
    """
    Curtailment cost should be zero if the zone is included in an RPS target
    because in this case the cost of curtailment will be calculated endogenously

    Args:
      model:
      zone:
      curtailment_cost:
      period:

    Returns:

    """
    if model.include_RPS:
        if resource in model.RPS_ELIGIBLE_RESOURCES_ALL:
            return curtailment_cost == 0

    return curtailment_cost >= 0

def emissions_groups_init(model):
    return sorted(set(
        val[2] for val in model.RESOURCE_EMISSIONS_GROUPS_MAP
    ))


def emissions_group_type_init(model, group):
    """
    Return the type of emissions for the current group, and validate that there is only one type of emissions for
    this group. This parameter is only used for validation and in reporting.
    """

    emission_type_list_cur_group = [
        emission_type
        for resource, emission_type, emissions_group in model.RESOURCE_EMISSIONS_GROUPS_MAP
        if emissions_group == group
    ]
    emission_types_cur_group = set(emission_type_list_cur_group)
    if len(emission_types_cur_group) > 1:
        raise ValueError("Each emissions group must cover only one emissions type")

    return list(emission_types_cur_group)[0]


def resource_emission_factor_denominator_validate(model, resource_emission_factor_denominator_mmbtu_or_mwh, resource, emission_type):
    """
    Resource emission factor denominator must be either mmbtu on mwh for thermal resources, or mwh for other resources
    """
    if (resource in model.THERMAL_RESOURCES) & (resource_emission_factor_denominator_mmbtu_or_mwh in ["mmbtu", "mwh"]):
        return True
    elif resource_emission_factor_denominator_mmbtu_or_mwh == "mwh":
        return True
    else:
        return False


def resource_emission_factor_validate_rule(model, group, fuel):
    """
    Loop through all resources of a given fuel for the current group, and check if there is a discrepancy in zero
    vs nonzero emission factors.

    The first for loop finds all the resource emission factors under the current group, and enters them in a dict,
    "resource_emission_factor_cur_fuel".

    The second if statement block first checks if any of the emission factors for the current group and fuel
    are nonzero (i.e., there are emissions defined for the current group). Then, it checks if any one of the emission
    factors is nonzero (i.e. at least one emission factor is zero when others are nonzero).
    """

    resource_emission_factors_cur_fuel = dict()

    for resource_cur in model.THERMAL_RESOURCES:
        if model.fuel[resource_cur] == fuel:
            emission_factor_current_resource_and_group = [
                model.resource_emission_factor[resource, emission_type]
                for resource, emission_type, emission_group in model.RESOURCE_EMISSIONS_GROUPS_MAP
                if (resource == resource_cur) & (emission_group == group)
            ]
            if emission_factor_current_resource_and_group == []:
                continue
            elif len(emission_factor_current_resource_and_group) > 1:
                raise ValueError("Resource {} has more than one mapping to group {}".format(resource_cur,group))
            else:
                resource_emission_factors_cur_fuel[resource_cur] = emission_factor_current_resource_and_group[0]

    resource_emission_factors_cur_fuel_list = list(resource_emission_factors_cur_fuel.values())

    if sum(resource_emission_factors_cur_fuel_list) > 0:
        if sum([x == 0 for x in resource_emission_factors_cur_fuel_list]) > 0:
            model.logger.warning(
                "Some resources with the fuel {} for policy {} have an emission factor of zero, while others are nonzero.\n"
                "The emission factors for this fuel under this policy are:\n{}"
                .format(fuel, group, resource_emission_factors_cur_fuel)
            )
    return True


def emission_factor_negative_direction_validate(model, emission_factor_negative_direction, line, group, period):
    """
    Validate that if a nonzero negative direction emission factor for a line-group-period combo is entered, the
    positive direction emission factor is zero. In other words, only one of these (positive and negative direction
    import rates) can be nonzero at a time.
    """
    if (
            emission_factor_negative_direction > 0
    ) & (
            model.emission_factor_positive_direction_per_mwh[line, group, period] > 0
    ):
        return False
    else:
        return True

# ### Reserves ### #


    # ############ Multi-day energy sharing params ###########


# energy sharing group set
def _energy_sharing_interval_init(model):
    """
    Get the set of values for interday_energy_sharing_interval_id (the set of energy sharing groups)
    """
    interday_energy_sharing_interval_ids = list()
    for day in model.DAYS:
        interday_energy_sharing_interval_ids.append(model.interday_energy_sharing_interval_id[day])

    return sorted(set(interday_energy_sharing_interval_ids))


def new_installed_capacity_tracking_rule(model, new_build_resource, period):
    """Track how much new capacity we have of each resource in each period.
    For each period, calculate how much capacity was installed in prior periods up to and including the current period.
    Used in results reporting.

    Args:
      model:
      period:
      new_build_resource:

    Returns:

    """
    new_installed_capacity = float()
    for vintage in model.VINTAGES:
        if vintage <= period:
            new_installed_capacity += model.Build_Capacity_MW[new_build_resource, vintage]

    return new_installed_capacity


def existing_tx_reliability_contribution_rule(model, reliability, tx_line, period):
    """
    min_flow_planned_mw is nonPositive; max_flow_planned_mw is nonNegative
    """
    if model.transmission_direction[reliability, tx_line, period] == 1:
        return (model.max_flow_planned_mw[tx_line] *
                model.existing_tx_reliability_contribution[reliability, tx_line, period])
    elif model.transmission_direction[reliability, tx_line, period] == -1:
        return (model.min_flow_planned_mw[tx_line] * -1 *
                model.existing_tx_reliability_contribution[reliability, tx_line, period])


def new_tx_reliability_contribution_rule(model, reliability, tx_line, period):
    """
    For new transmission, there is no direction to the tx line
    """
    if tx_line in model.TRANSMISSION_LINES_NEW:
        return (model.New_Tx_Total_Installed_Capacity_MW[tx_line, period] *
                model.new_tx_reliability_contribution[reliability, tx_line, period])
    else:
        return 0


def cumulative_new_storage_charge_capacity_build_tracking_rule(model, new_build_storage_resource, period):
    """
    Track how much new storage charge capacity we have of each resource in each period.
    For each period, calculate how much storage charge capacity was built in prior periods up to and including the
    current period.
    Used in calculating operational charging capacity.
    """
    total_charge_capacity_build = float()
    for vintage in model.VINTAGES:
        if vintage <= period:
            total_charge_capacity_build += model.New_Storage_Charge_Capacity_Build_MW[new_build_storage_resource,
                                                                                      vintage]

    return total_charge_capacity_build


def retire_new_capacity_period_total_rule(model, resource, period):
    """Track how much new capacity has been retired up to and including the current period
    Used in results reporting.

    Args:
      model:
      resource:
      period:

    Returns:

    """
    retire_new_capacity = float()
    for vintage in model.VINTAGES:
        if period >= vintage:
            retire_new_capacity += model.Retire_New_Capacity_By_Vintage_Cumulative_MW[resource, period, vintage]

    return retire_new_capacity


def operational_new_capacity_in_period_by_vintage_rule(model, resource, period, vintage):
    """The amount of new capacity that is still in operation (not retired) in a period
    following installation in a previous vintage. Elsewhere retirements when period = vintage are prohibited.

    Args:
      model:
      resource:
      period:
      vintage:

    Returns:

    """

    if resource in model.CAN_RETIRE_RESOURCES_NEW:
        retired_capacity = model.Retire_New_Capacity_By_Vintage_Cumulative_MW[resource, period, vintage]
    else:
        retired_capacity = 0

    return model.Build_Capacity_MW[resource, vintage] - retired_capacity


def operational_new_capacity_in_period_rule(model, resource, period):
    """Track how much new capacity is still operational (has not been retired) in each period.

    Args:
      model:
      resource:
      period:

    Returns:

    """
    operational_new_capacity = float()
    for vintage in model.VINTAGES:
        if period >= vintage:
            operational_new_capacity += \
                model.Operational_New_Capacity_In_Period_By_Vintage_MW[resource, period, vintage]

    return operational_new_capacity


def cumulative_retired_planned_capacity_rule(model, resource, period):
    """Track how much planned capacity has been retired up to and including the current period
    Used in results reporting.

    Args:
      model:
      period:
      resource:

    Returns:

    """
    cumulative_retired_planned_capacity = float()
    for p in model.PERIODS:
        if p <= period:
            cumulative_retired_planned_capacity += model.Retire_Planned_Capacity_MW[resource, p]

    return cumulative_retired_planned_capacity


# Define the planned capacity actually operational in each period
# This will only diverge from the planned capacity input values if retirements are allowed
def operational_planned_capacity_in_period_rule(model, resource, period):
    """

    Args:
      model:
      resource:
      period:

    Returns:

    """

    if resource in model.CAN_RETIRE_RESOURCES:
        cumulative_planned_retirements_mw = model.Retire_Planned_Capacity_Cumulative_MW[resource, period]
    else:
        # if endogenous retirements aren't allowed, the planned capacity can still be reduced over time
        # via the schedule of planned capacity subtractions
        cumulative_planned_retirements_mw = model.cumulative_planned_capacity_subtractions_mw[resource, period]

    # planned capacity in the first period is an "addition"
    return model.cumulative_planned_capacity_additions_mw[resource, period] - cumulative_planned_retirements_mw


def operational_capacity_tracking_rule(model, resource, period):
    """Track how much capacity of each resource is operational in each period,
    including planned capacity, new builds, and retirements

    Args:
      model:
      period:
      resource:

    Returns:

    """
    if resource in model.NEW_BUILD_RESOURCES:
        operational_new_capacity = model.Operational_New_Capacity_MW[resource, period]
    else:
        operational_new_capacity = 0

    return model.Operational_Planned_Capacity_MW[resource, period] + operational_new_capacity


def available_capacity_def(model, resource, timepoint):
    """Define the amount of capacity available for dispatch in a given timepoint,
    which is the Operational Capacity in the period, derated by outages.

    Args:
      model:
      resource:
      timepoint:

    Returns:

    """

    return model.Operational_Capacity_MW[resource, model.period[timepoint]] \
           * model.available_capacity_fraction[resource, timepoint]


def max_increase_load_in_timepoint_rule(model, resource, timepoint):
    """
    Define the amount of max increase load for dispatch in a given timepoint.
    New build resources include the new capacity in the calculation.
    """

    if resource in model.RESOURCES_WITH_ENERGY_IN_STORAGE \
            and model.planned_storage_charge_capacity_mw[resource, model.period[timepoint]] != float("inf"):
        model.logger.warning(
            "The parameter `planned_storage_charge_capacity_mw` is deprecated. "
            "Use `max_increase_load_fraction` if there is a desire to limit storage charging capacity."
        )
        return model.Operational_Charging_Capacity_MW[resource, model.period[timepoint]] \
               * model.available_capacity_fraction[resource, timepoint] \
               * model.max_increase_load_fraction[resource, timepoint]
    else:
        return model.Available_Capacity_In_Timepoint_MW[resource, timepoint] \
                   * model.max_increase_load_fraction[resource, timepoint]


def min_increase_load_in_timepoint_rule(model, resource, timepoint):
    """
    Define a minimum load in a given timepoint.
    New build resources include the new capacity in the calculation.
    """

    if resource in model.RESOURCES_WITH_ENERGY_IN_STORAGE \
            and model.planned_storage_charge_capacity_mw[resource, model.period[timepoint]] != float("inf"):
        model.logger.warning(
            "The parameter `planned_storage_charge_capacity_mw` is deprecated. "
            "Use `max_increase_load_fraction` if there is a desire to limit storage charging capacity."
        )
        return model.Operational_Charging_Capacity_MW[resource, model.period[timepoint]] \
               * model.available_capacity_fraction[resource, timepoint] \
               * model.min_increase_load_fraction[resource, timepoint]
    else:
        return model.Available_Capacity_In_Timepoint_MW[resource, timepoint] \
                   * model.min_increase_load_fraction[resource, timepoint]


def operational_charging_capacity_rule(model, resource, period):
    """
    Define operational charging capacity. Equal to operational capacity if not specifying different size charge from
    discharge. Else, set equal to the total available charge capacity, including new build if applicable.

    While Cumulative_New_Storage_Charge_Capacity_Build_MW will exist in all cases, it will only show up in the
    optimization if there is a planned storage charge capacity defined (which could equal zero)
    """

    if model.planned_storage_charge_capacity_mw[resource, period] == float("inf"):
        return model.Operational_Capacity_MW[resource, period]
    else:
        if resource not in model.NEW_BUILD_RESOURCES:
            return model.planned_storage_charge_capacity_mw[resource, period]
        else:
            return (model.planned_storage_charge_capacity_mw[resource, period] +
                    model.Cumulative_New_Storage_Charge_Capacity_Build_MW[resource, period])


# Storage energy capacity expressions
def new_storage_installed_energy_capacity_tracking_rule(model, new_build_storage_resource, current_period):
    """Track how much new energy capacity we have of each storage resource in each period.
    For each period, calculate how much energy capacity was installed in prior periods up to and including the current
    period.

    Args:
      model:
      current_period:
      new_build_storage_resource:

    Returns:

    """
    new_storage_installed_energy_capacity = float()
    for vintage in model.VINTAGES:
        if vintage <= current_period:
            new_storage_installed_energy_capacity += \
                model.Build_Storage_Energy_Capacity_MWh[new_build_storage_resource, vintage]

    return new_storage_installed_energy_capacity


def total_storage_installed_energy_capacity_tracking_rule(model, storage_resource, period):
    """Track how much total energy capacity we have of each storage resource in each period,
    including planned capacity and new builds.
    For each period, calculate the sum of planned capacity in the period and how much capacity was installed in prior
    periods up to and including the current period.

    Args:
      model:
      period:
      storage_resource:

    Returns:

    """
    if storage_resource in model.NEW_BUILD_RESOURCES:
        cumulative_new_build_mwh = model.Cumulative_New_Storage_Energy_Capacity_MWh[storage_resource, period]
    else:
        cumulative_new_build_mwh = 0

    return model.planned_storage_energy_capacity_mwh[storage_resource, period] + cumulative_new_build_mwh


def reliability_capacity_delivered_mw_rule(model, resource, period):
    """Returns reliability capacity as either A) Deliverability_Capacity_MW
    or B) Operational_Capacity_MW
    Used in reliability_preNQC_power_rule, reliability_duration_MWh_rule and
    reliability_preNQC_conventional_dr_annual_rule

    Args:
      model:
      deliverability_zone:
      period:

    Returns:

    """
    if model.deliverability and resource in model.DELIVERABILITY_RESOURCES:
        return model.Delivered_Capacity_MW[resource, period]
    else:
        return model.Operational_Capacity_MW[resource, period]



def prestart_to_start_definition(model, resource, timepoint):
    """A unit starts (adds to Commit_Units) min_down_time_hours after it was pre-started

    Args:
      model:
      timepoint:
      resource:

    Returns:

    """
    min_down_time = model.min_down_time_hours[resource]

    # Resolve days are circular, so if adding min_down_time moves past the end of the day, wrap around to the start
    if timepoint + min_down_time > model.last_timepoint_of_day[model.period[timepoint], model.day[timepoint]]:
        day_wrap = model.timepoints_per_day
    else:
        day_wrap = 0

    return model.Start_Units[resource, timepoint + min_down_time - day_wrap]


def preshut_down_to_shut_down_definition(model, resource, timepoint):
    """A unit shuts down (subtracts from to Commit_Units) min_up_time_hours after it was preshut down

    Args:
      model:
      timepoint:
      resource:

    Returns:

    """
    min_up_time = model.min_up_time_hours[resource]

    # Resolve days are circular, so if adding min_down_time moves past the end of the day, wrap around to the start
    if timepoint + min_up_time > model.last_timepoint_of_day[model.period[timepoint], model.day[timepoint]]:
        day_wrap = model.timepoints_per_day
    else:
        day_wrap = 0

    return model.Shut_Down_Units[resource, timepoint + min_up_time - day_wrap]


def define_starting_units(model, resource, timepoint):
    """Track how many units are in the process of starting in each timepoint
    by summing over previous timepoints within the min down time for each resource.
    This variable will then be used to determine commitment from timepoint to timepoint
    (in conjunction with the variable tracking units scheduled to shut down).

    Args:
      model:
      timepoint:
      resource:

    Returns:

    """
    starting_not_yet_committed_units = float()

    min_down_time = model.min_down_time_hours[resource]

    # Days are treated as circular, so if min down time puts us across the boundary (hour 24 to hour 1),
    # add 24 hours to avoid negative numbers (for example if we are in timepoint 3 and have a 5-hour min down time,
    # we'll need to look at units scheduled to start in the four previous hours,
    # namely hours (3+24)-4=23, (3+24)-3=24, 1, and 2.
    if timepoint - min_down_time < model.first_timepoint_of_day[model.period[timepoint], model.day[timepoint]]:
        previous_timepoints = list()
        for pt in range(1, min_down_time):
            if timepoint - pt >= model.first_timepoint_of_day[model.period[timepoint], model.day[timepoint]]:
                previous_timepoints.append(timepoint - pt)
            else:
                previous_timepoints.append(timepoint + model.timepoints_per_day - pt)

        for t in previous_timepoints:
            starting_not_yet_committed_units += model.PreStart_Units[resource, t]
    else:
        for t in range(1, min_down_time):
            starting_not_yet_committed_units += model.PreStart_Units[resource, timepoint - t]

    return starting_not_yet_committed_units


def define_shutting_down_units(model, resource, timepoint):
    """Track how many units are in the process of shutting down in each timepoint
    by summing over previous timepoints within the min up time for each resource.
    This variable will then be used to determine commitment from timepoint to timepoint
    (in conjunction with the variable tracking units scheduled to start up).

    Args:
      model:
      timepoint:
      resource:

    Returns:

    """

    shutting_down_still_committed_units = float()

    min_up_time = model.min_up_time_hours[resource]

    # We'll treat the day as circular, so if min up time puts us across the boundary (hour 24 to hour 1), we'll have
    # to add 24 hours to avoid negative numbers (for example if we are in timepoint 3 and have a 5-hour min up time,
    # we'll need to look at units scheduled to shut down in the four previous hours,
    # namely hours (3+24)-4=23, (3+24)-3=24, 1, and 2.
    if timepoint - min_up_time < model.first_timepoint_of_day[model.period[timepoint], model.day[timepoint]]:
        previous_timepoints = list()
        for pt in range(1, min_up_time):
            if timepoint - pt >= model.first_timepoint_of_day[model.period[timepoint], model.day[timepoint]]:
                previous_timepoints.append(timepoint - pt)
            else:
                previous_timepoints.append(timepoint + model.timepoints_per_day - pt)

        for t in previous_timepoints:
            shutting_down_still_committed_units += model.PreShut_Down_Units[resource, t]
    else:
        for t in range(1, min_up_time):
            shutting_down_still_committed_units += model.PreShut_Down_Units[resource, timepoint - t]

    return shutting_down_still_committed_units


def fully_operational_units_def_rule(model, resource, timepoint):
    """The number of units that are fully operational is the number that aren't in their start or stop sequence,
    which is currently to hold the unit at PMin for one timestep.
    Shut_Down_Units from the next timepoint is used because this variable is non-zero in
    the timestep *after* the unit has finished its stop sequence
    Fully_Operational_Units is only defined for the set UNIT_COMMITMENT_RAMP_LIMITED_RESOURCES because
    if a unit can ramp quickly then all committed units are fully operational

    Args:
      model:
      timepoint:
      resource:

    Returns:

    """
    return model.Commit_Units[resource, timepoint] \
           - model.Start_Units[resource, timepoint] \
           - model.Shut_Down_Units[resource, model.next_timepoint[timepoint]]


def define_start_capacity(model, resource, timepoint):
    """Calculate MW of capacity started in timepoint

    Args:
      model:
      resource: unit commitment resource
      timepoint:

    Returns:

    """
    return model.Start_Units[resource, timepoint] * model.unit_size_mw[resource]


def define_commit_capacity(model, resource, timepoint):
    """Calculate MW of capacity committed in timepoint

    Args:
      model:
      resource: unit commitment resource
      timepoint:

    Returns:

    """
    return model.Commit_Units[resource, timepoint] * model.unit_size_mw[resource]


def define_sales_after_ee_mwh(model, policy, period):
    """Reduce the RPS target when the model builds EE
    The capacity of EE resources is in average MW - multiplying by hours_per_year gives yearly MWh

    Args:
      model:
      period:

    Returns:

    """
    total_ee_mwh = float()
    for ee_program in model.EE_PROGRAMS:
        if ee_program in model.RPS_ELIGIBLE_RESOURCES[policy]:
            total_ee_mwh += (model.Operational_Capacity_MW[ee_program, period] *
                             model.hours_per_year *
                             model.rps_eligible_fraction[policy, ee_program])

    return model.retail_sales_mwh[policy, period] - total_ee_mwh


def define_rps_target(model, policy, period):
    """Reduce the RPS target when the model builds EE
    The capacity of EE resources is in average MW - multiplying by hours_per_year gives yearly MWh

    Args:
      model:
      period:

    Returns:

    """
    return model.Retail_Sales_After_EE_MWh[policy, period] * model.rps_fraction_of_retail_sales[policy, period]


def define_ee_program_load_in_timepoint(model, resource, timepoint):
    """Derive EE program load reduction contribution based on the total EE program investment and EE impact shape
    BTM = behind the meter - there are avoided T & D losses from EE investments.

    Args:
      model:
      timepoint:
      resource:

    Returns:

    """
    return model.Operational_Capacity_MW[resource, model.period[timepoint]] \
           * model.ee_shapes_btm_mwh_per_amw[resource, model.day[timepoint], model.hour_of_day[timepoint]] \
           * (1 + model.ee_t_and_d_losses_fraction[resource])


def downward_lf_subhourly_energy_rule(model, zone, timepoint):
    """Calculate energy dispatch by variable resources providing LF down.
    # Allocate the subhourly energy to each constituting zone of the LOAD_Following zones, according to input load

    Args:
      model:
      timepoint:

    Returns:

    """
    # Tally up the total load in all LF zones
    total_load_in_LF_zone = 0
    for LF_zone in model.LOAD_FOLLOWING_ZONES:
        total_load_in_LF_zone += model.input_load_mw[LF_zone, timepoint]

    # calculate the portion of load taken up by the current zone
    lf_zone_load_fraction = model.input_load_mw[zone, timepoint]/total_load_in_LF_zone

    # Prorate the subhourly energy expenditure based on input_load
    return (model.lf_reserve_dispatch_fraction * model.Variable_Resource_Provide_Downward_LF_MW[timepoint]
            *lf_zone_load_fraction)


def upward_lf_subhourly_energy_rule(model, zone, timepoint):
    """Calculate energy dispatch by variable resources providing LF up.

    Args:
      model:
      timepoint:

    Returns:

    """
    # Tally up the total load in all LF zones
    total_load_in_LF_zone = 0
    for LF_zone in model.LOAD_FOLLOWING_ZONES:
        total_load_in_LF_zone += model.input_load_mw[LF_zone, timepoint]

    # calculate the portion of load taken up by the current zone
    lf_zone_load_fraction = model.input_load_mw[zone, timepoint] / total_load_in_LF_zone

    # Prorate the subhourly energy expenditure based on input_load
    return (model.lf_reserve_dispatch_fraction * model.Variable_Resource_Provide_Upward_LF_MW[timepoint]
            * lf_zone_load_fraction)


def define_transmit_power(model, transmission_line, timepoint):
    """Expression that adds unspecified flows along each transmission line along to specified imports and exports
    to arrive at the net flow along the line.
    This represents the metered flow at the boundary between balancing areas.

    Args:
      model:
      timepoint:
      transmission_line:

    Returns:

    """
    transmission_reserved_for_dedicated_resources = 0.0

    if model.resource_use_tx_capacity:
        for resource_line_pair in model.RESOURCE_TX_IDS:
            resource = model.dedicated_import_resource[resource_line_pair]
            if transmission_line == model.tx_line_used[resource_line_pair]:
                if resource in model.RESOURCES_THAT_INCREASE_LOAD:
                    charge_storage = model.Increase_Load_MW[resource, timepoint]
                else:
                    charge_storage = 0.0

                transmission_reserved_for_dedicated_resources += \
                    (model.Provide_Power_MW[resource, timepoint] - charge_storage) \
                    * model.resource_tx_direction[resource_line_pair]

    return model.Transmit_Power_Unspecified_MW[transmission_line, timepoint] \
           + transmission_reserved_for_dedicated_resources


# Expressions for hurdle rate costs associated with unspecified imports in each timepoint.
def define_hurdle_cost_per_timepoint_positive_direction(model, line, timepoint):
    """

    Args:
      model:
      line:
      timepoint:

    Returns:

    """

    return model.Transmit_Power_Unspecified_Positive_Direction_MW[line, timepoint] \
           * model.positive_direction_hurdle_rate_per_mw[line, model.period[timepoint]]


def define_hurdle_cost_per_timepoint_negative_direction(model, line, timepoint):
    """

    Args:
      model:
      line:
      timepoint:

    Returns:

    """

    return model.Transmit_Power_Unspecified_Negative_Direction_MW[line, timepoint] \
           * model.negative_direction_hurdle_rate_per_mw[line, model.period[timepoint]]


# Expressions for emissions attributed to unspecified imports in each timepoint. Cannot be negative.
# Transmit_Power_Unspecified is used here because dedicated imports
# are modeled as supplying power to the "to" zone directly
def define_unspecified_import_emissions_per_timepoint_positive_direction(model, line, group, timepoint):
    """

    Args:
      model:
      line:
      group:
      timepoint:

    Returns:

    """

    return model.Transmit_Power_Unspecified_Positive_Direction_MW[line, timepoint] \
        * model.emission_factor_positive_direction_per_mwh[line, group, model.period[timepoint]]


def define_unspecified_import_emissions_per_timepoint_negative_direction(model, line, group, timepoint):
    """

    Args:
      model:
      line:
      group:
      timepoint:

    Returns:

    """

    return model.Transmit_Power_Unspecified_Negative_Direction_MW[line, timepoint] \
        * model.emission_factor_negative_direction_per_mwh[line, group, model.period[timepoint]]


def define_unspecified_import_emissions(model, group, period):
    """Expression that sums the emissions from imports into the emissions target area for each period

    Args:
      model:
      period:

    Returns:

    """
    unspecified_import_units_per_year = float()

    for timepoint in model.TIMEPOINTS:
        if model.period[timepoint] == period:
            for line in model.TRANSMISSION_LINES:
                unspecified_import_units_per_year += \
                    (model.Unspecified_Imports_Units_Per_Timepoint_Positive_Direction[line, group, timepoint]
                    + model.Unspecified_Imports_Units_Per_Timepoint_Negative_Direction[line, group, timepoint]) \
                * model.day_weight[model.day[timepoint]]

    return unspecified_import_units_per_year


def sum_zonal_increase_load_rule(model, zone, timepoint):
    """
    sum Increase_Load_MW for all resources in a zone
    :param model:
    :param zone:
    :param timepoint:
    :return:
    """

    increase_load = float()
    for resource in model.RESOURCES_THAT_INCREASE_LOAD:
        if model.zone[resource] == zone:
            increase_load += model.Increase_Load_MW[resource, timepoint]

    return increase_load


def sum_zonal_provide_power_rule(model, zone, timepoint):
    """
    sum Provide_Power_MW for all resources in a zone
    :param model:
    :param zone:
    :param timepoint:
    :return:
    """

    provide_power = float()
    for resource in model.RESOURCES_WITH_PROVIDE_POWER:
        if model.zone[resource] == zone:
            provide_power += model.Provide_Power_MW[resource, timepoint]

    return provide_power


def sum_zonal_net_provide_power_rule(model, zone, timepoint):
    """
    get the net provide power in each zone
    :param model:
    :param zone:
    :param timepoint:
    :return:
    """

    return model.Zonal_Provide_Power_MW[zone, timepoint] - model.Zonal_Increase_Load_MW[zone, timepoint]


def define_blend_fuel_generation_mw(model, resource, timepoint):
    """This expression defines how much power is generated from blend fuel for the purpose of RPS accounting.
    This rule is effective during most conditions (not only Pmax), although in terms of limitation:
    1. If the unit went through a lot of start stop, its heat rate would deviate from average heat rate.
    2. In rare cases when blend_fuel provided almost all fuel consumption for a year, we might see blend_fuel_gen
    exceeding provide_power by a small (<1%) amount.

    Args:
      model:
      timepoint:
      resource:

    Returns:

    """

    return model.Blend_Fuel_Consumption_MMBtu[resource, timepoint] / model.avg_heat_rate_mmbtu_per_mwh[resource]


def define_yearly_fuel_consumption(model, fuel, period):
    """Expression that represents the yearly consumption of fuel.

    Fuels can be burnt as either "base" fuel or the "blend" fuel, and we want to capture both here.

    The naming is a bit confusing, since here we are summing up to the FUEL, whereas elsewhere,
    we also sum up to the fuel burned by each RESOURCE.

    Args:
      model:
      fuel:
      period:

    Returns:

    """
    yearly_base_fuel_consumption = 0.0
    yearly_blend_fuel_consumption = 0.0
    for timepoint in model.TIMEPOINTS:
        if model.period[timepoint] == period:
            for resource in model.THERMAL_RESOURCES:
                # Add up base fuel consumption
                if model.fuel[resource] == fuel:
                    yearly_base_fuel_consumption  += (
                        model.Base_Fuel_Consumption_MMBtu[resource, timepoint] *
                        model.day_weight[model.day[timepoint]]
                    )
                # Add up blend fuel consumption
                if resource in model.BLEND_FUEL_RESOURCES:
                    if model.blend_fuel[model.fuel[resource]] == fuel:
                        yearly_blend_fuel_consumption += (
                            model.Blend_Fuel_Consumption_MMBtu[resource, timepoint] *
                            model.day_weight[model.day[timepoint]]
                        )

    return yearly_base_fuel_consumption + yearly_blend_fuel_consumption


def define_yearly_blend_fuel_rps_generation(model, policy, period):
    """derived variable that represents the yearly generation of electricity from blend fuel
    that counts towards the RPS target

    Args:
      model:
      policy:
      period:

    Returns:

    """
    yearly_blend_fuel_RPS_generation = 0.0
    for resource in (model.BLEND_FUEL_RESOURCES & model.RPS_ELIGIBLE_RESOURCES[policy]):
        for timepoint in model.TIMEPOINTS:
            if model.period[timepoint] == period:
                yearly_blend_fuel_RPS_generation += (model.Blend_Fuel_Generation_MW[resource, timepoint]
                                                     * model.day_weight[model.day[timepoint]]
                                                     * model.rps_eligible_fraction[policy, resource])

    return yearly_blend_fuel_RPS_generation


def define_start_fuel_in_timepoint(model, resource, timepoint):
    """Calculates start fuel consumption for each unit commitment resource in each timepoint.

    Args:
      resource: unit commitment resource
      timepoint:
      model:

    Returns:

    """
    return model.Start_Capacity_MW[resource, timepoint] * model.start_fuel_mmbtu_per_mw[resource]


def define_base_fuel_cost(model, resource, timepoint):
    """Cost of fuel for each thermal resource for each timepoint. Base fuel

    Args:
      model:
      resource:
      timepoint:

    Returns:

    """
    return (model.Base_Fuel_Consumption_MMBtu[resource, timepoint] *
            model.fuel_price_per_mmbtu[model.fuel[resource],
                                       model.period[timepoint], model.month[timepoint]])


def define_blend_fuel_cost(model, resource, timepoint):
    """Cost for fuel for each thermal resource for each timepoint. Blend Fuel

    Args:
      model:
      resource:
      timepoint:

    Returns:

    """
    return (model.Blend_Fuel_Consumption_MMBtu[resource, timepoint]*
            model.fuel_price_per_mmbtu[model.blend_fuel[model.fuel[resource]],
                                       model.period[timepoint], model.month[timepoint]])


def define_total_fuel_cost(model, resource, timepoint):
    """Cost for fuel for each thermal resource for each timepoint. Summed with base fuel and blend fuel

    Args:
      model:
      resource:
      timepoint:

    Returns:

    """
    if resource in model.BLEND_FUEL_RESOURCES:
        blend_fuel_costs = model.Blend_Fuel_Cost[resource,timepoint]
    else:
        blend_fuel_costs = 0.0

    return model.Base_Fuel_Cost[resource,timepoint] + blend_fuel_costs

def define_resource_timepoint_emissions(model, resource, emission_type, timepoint):
    """Hourly emissions for each resource in each timepoint
    The variables Base_Fuel_Consumption_MMBtu and Blend_Fuel_Consumption_MMBtu include fuel consumption from normal operations as well as starting and reserve dispatch.

    Args:
      model:
      resource:
      group:
      timepoint:

    Returns:

    """
    if model.resource_emission_factor_denominator_mmbtu_or_mwh[resource, emission_type] == "mmbtu":

        if resource in model.THERMAL_RESOURCES:
            base_fuel_emissions = (model.Base_Fuel_Consumption_MMBtu[resource, timepoint]
                                    * model.resource_emission_factor[resource, emission_type])
        else:
            base_fuel_emissions = 0

        if resource in model.BLEND_FUEL_RESOURCES:
            blend_fuel_emissions = (model.Blend_Fuel_Consumption_MMBtu[resource, timepoint]
                                    * model.resource_emission_factor_blended_fuel[resource, emission_type])
        else:
            blend_fuel_emissions = 0.0

        resource_emissions_in_timepoint = base_fuel_emissions + blend_fuel_emissions

    elif model.resource_emission_factor_denominator_mmbtu_or_mwh[resource, emission_type] == "mwh":

        power_output_mwh = model.Provide_Power_MW[resource, timepoint]

        # Regulation reserve dispatch
        if resource in model.REGULATION_RESERVE_RESOURCES:
            power_output_mwh += (model.Provide_Upward_Reg_Net_MW[resource, timepoint]
                                 - model.Provide_Downward_Reg_Net_MW[resource, timepoint]) \
                                 * model.reg_dispatch_fraction

        # Load-following reserve dispatch
        if resource in model.LOAD_FOLLOWING_RESERVE_RESOURCES:
            power_output_mwh += (model.Provide_LF_Upward_Reserve_Net_MW[resource, timepoint]
                                 - model.Provide_LF_Downward_Reserve_Net_MW[resource, timepoint]) \
                                 * model.lf_reserve_dispatch_fraction

        resource_emissions_in_timepoint = power_output_mwh * model.resource_emission_factor[resource, emission_type]
    else:
        raise ValueError("emission factor denominator units must be mmbtu or mwh")

    return resource_emissions_in_timepoint


def define_resource_emissions(model, group, period):
    """Expression that in each period that sums the emissions from resources covered by the current group

    Args:
      model:
      group:
      period:

    Returns:

    """
    return sum(
                sum(
                    model.Resource_Emissions_Units_Per_Timepoint[resource, emission_type, timepoint]
                    * model.day_weight[model.day[timepoint]]
                    for resource, emission_type, emissions_group in model.RESOURCE_EMISSIONS_GROUPS_MAP
                    if emissions_group == group
                )
                for timepoint in model.TIMEPOINTS
                if model.period[timepoint] == period
            )


##################################
# ##### OBJECTIVE FUNCTION ##### #
##################################
def define_variable_cost_in_timepoint(model, resource, timepoint):
    """Derived variable that represents the variable cost of each resource in each timepoint (not weighted or discounted)
    All resources are included here, even if variable cost is zero

    Args:
      model:
      timepoint:
      resource:

    Returns:

    """

    # variable cost to increase load
    if resource in model.RESOURCES_THAT_INCREASE_LOAD:
        load_cost = model.Increase_Load_MW[resource, timepoint] \
                    * model.variable_cost_increase_load_per_mwh[resource, timepoint]
    else:
        load_cost = 0.0
    # variable cost to provide power
    # define the variable cost of EE programs behind the meter.
    if resource in model.EE_PROGRAMS and model.allow_ee_investment:
        power_var = model.Operational_Capacity_MW[resource, model.period[timepoint]] \
                    * model.ee_shapes_btm_mwh_per_amw[resource, model.day[timepoint], model.hour_of_day[timepoint]]
    elif resource in model.RESOURCES_WITH_PROVIDE_POWER:
        power_var = model.Provide_Power_MW[resource, timepoint]
    else:
        power_var = 0.0
    power_cost = power_var * model.variable_cost_per_mwh[resource, timepoint]

    return load_cost + power_cost


def define_start_cost_in_timepoint(model, resource, timepoint):
    """Calculates start costs for each unit commitment resource in each timepoint.

    Args:
      resource: unit commitment resource
      timepoint:
      model:

    Returns:

    """
    return model.Start_Capacity_MW[resource, timepoint] * model.startup_cost_per_mw[resource]


def define_shutdown_cost_in_timepoint(model, resource, timepoint):
    """Calculates shutdown costs for each unit commitment resource in each timepoint.

    Args:
      resource: unit commitment resource
      timepoint:
      model:

    Returns:

    """
    return model.Shut_Down_Units[resource, timepoint] \
           * model.unit_size_mw[resource] \
           * model.shutdown_cost_per_mw[resource]


def capital_cost_annual_in_period_dollars_rule(model, resource, period):
    """

    Args:
      model:
      resource:
      period:

    Returns:

    """
    # add levelized capital costs for all resource additions from previous periods and the current period
    # Note that capital cost inputs are in terms of kW and kWh instead of MW or MWh
    # because industry practice is to express levelized capital costs as $/kW-yr

    capital_cost_total = 0.0

    for vintage in model.VINTAGES:
        if period >= vintage:
            capital_cost_total += model.Build_Capacity_MW[resource, vintage] \
                                  * model.capital_cost_per_kw_yr[resource, vintage] * 10 ** 3

            # if storage, also add energy capacity costs
            if resource in model.NEW_BUILD_STORAGE_RESOURCES:
                capital_cost_total += model.Build_Storage_Energy_Capacity_MWh[resource, vintage] \
                                      * model.energy_storage_cost_per_kwh_yr[resource, vintage] * 10 ** 3 \
                                      + model.New_Storage_Charge_Capacity_Build_MW[resource,vintage]\
                                      * model.energy_storage_charge_capacity_cost_per_kw_yr[resource, vintage] * 10 ** 3

    return capital_cost_total


def fixed_o_and_m_annual_in_period_dollars_rule(model, resource, period):
    """

    Args:
      model:
      resource:
      period:

    Returns:

    """

    fixed_om_total = 0.0

    # Fixed O&M for all planned capacity that is still operational.
    # For storage, this is only the power-related portion of Fixed O&M
    fixed_om_total += model.Operational_Planned_Capacity_MW[resource, period] \
                      * model.planned_capacity_fixed_o_and_m_dollars_per_kw_yr[resource, period] * 10 ** 3

    # Energy/duration-related portion of Fixed O&M for storage
    # Note: storage retirements are not currently allowed, so for storage energy capacity, planned = operational
    if resource in model.STORAGE_RESOURCES:
        fixed_om_total += model.planned_storage_energy_capacity_mwh[resource, period] \
            * model.planned_storage_energy_capacity_fixed_o_and_m_dollars_per_kwh_yr[resource, period]  * 10 ** 3

    # Fixed O&M for new build resources, which varies by installation date (vintage)
    if resource in model.NEW_BUILD_RESOURCES:
        for vintage in model.VINTAGES:
            if period >= vintage:
                fixed_om_total += model.Operational_New_Capacity_In_Period_By_Vintage_MW[resource, period, vintage] \
                                  * model.new_capacity_fixed_o_and_m_dollars_per_kw_yr[resource, vintage] * 10 ** 3

                # if storage, also add fixed O&M costs for maintaining energy capacity
                # Note: storage retirements are not currently allowed, so new = operational.
                if resource in model.NEW_BUILD_STORAGE_RESOURCES:
                    fixed_om_total += model.Build_Storage_Energy_Capacity_MWh[resource, vintage] \
                                      * model.new_energy_capacity_fixed_o_and_m_dollars_per_kwh_yr[resource, vintage] \
                                      * 10 ** 3

    return fixed_om_total


def total_cost_rule(model):
    """Total incremental cost of the new resources + system operational costs including variable and fuel costs,
    and cost of unserved energy and unserved reserves.

    Args:
      model:

    Returns:

    """

    # ### Annual (periodic) costs ### #

    # Capital and Fixed costs of existing and new resources
    capital_and_fixed_resource_costs = float()
    deliverability_cost = float()
    for period in model.PERIODS:
        for resource in model.RESOURCES:
            capital_and_fixed_resource_costs += \
                model.Fixed_OM_Annual_In_Period_Dollars[resource, period] * model.discount_factor[period]

            # for new resources, add capital costs
            if resource in model.NEW_BUILD_RESOURCES:
                capital_and_fixed_resource_costs += \
                    model.Capital_Cost_Annual_In_Period_Dollars[resource, period] * model.discount_factor[period]

        # Deliverability costs
        if model.deliverability:
            for deliverability_zone in model.DELIVERABILITY_ZONES:
                deliverability_cost += model.Deliverability_Cost_Annual_In_Period_Dollars[deliverability_zone, period] \
                    * model.discount_factor[period]

    if model.allow_tx_build:
        transmission_build_costs = float()
        for (p, v) in model.PERIOD_VINTAGES:
            for l in model.TRANSMISSION_LINES_NEW:
                transmission_build_costs += (model.New_Tx_Period_Vintage_Build_Cost[l, p, v] * model.discount_factor[p])

    # ### Hourly costs ### #
    variable_costs = float()
    fuel_costs = float()
    unit_start_costs = float()
    unit_shutdown_costs = float()
    unserved_energy_costs = float()
    overgen_costs = float()
    scheduled_curtailment_costs = float()
    upward_reg_violation_costs = float()
    downward_reg_violation_costs = float()
    upward_lf_reserve_violation_costs = float()
    downward_lf_reserve_violation_costs = float()
    spin_violation_costs = float()
    hurdle_rate_costs = float()

    for tmp in model.TIMEPOINTS:

        timepoint_weight = model.day_weight[model.day[tmp]] * model.discount_factor[model.period[tmp]]

        # ### Resource hourly costs ### #
        for r in model.RESOURCES:

            # Variable costs, excluding fuel
            variable_costs += model.Variable_Cost_In_Timepoint[r, tmp] * timepoint_weight

            # Thermal fuel cost
            if r in model.THERMAL_RESOURCES:
                fuel_costs += model.Total_Fuel_Cost[r, tmp] * timepoint_weight

            if r in model.UNIT_COMMITMENT_RESOURCES:
                unit_start_costs += (model.Start_Cost_In_Timepoint[r, tmp] * timepoint_weight)

                unit_shutdown_costs += (model.Shutdown_Cost_In_Timepoint[r, tmp] * timepoint_weight)

            # Scheduled curtailment costs
            # the assumed net cost of procuring additional renewable energy certificates
            # for each curtailable resource that isn't covered under an RPS target
            if r in model.CURTAILABLE_VARIABLE_RESOURCES:
                scheduled_curtailment_costs += model.curtailment_cost_per_mwh[r, model.period[tmp]] \
                    * model.Scheduled_Curtailment_MW[r, tmp]* timepoint_weight

        # ### Hurdle rates ### #
        # Only for unspecified power because dedicated imports and exports aren't charged hurdle rates
        for line in model.TRANSMISSION_LINES:
            hurdle_rate_costs += (model.Hurdle_Cost_Per_Timepoint_Positive_Direction[line, tmp]
                                  + model.Hurdle_Cost_Per_Timepoint_Negative_Direction[line, tmp]) \
                * timepoint_weight

        # add in hurdle rates for semi storage zones transmitted power
        if model.allow_semi_storage_zones:
            for ssz in model.SEMI_STORAGE_ZONES:
                hurdle_rate_costs += (model.SSZ_Positive_Transmit_Power_MW[ssz, tmp]
                                     * model.ssz_positive_direction_hurdle_rate_per_mw[ssz, model.period[tmp]]
                                     + model.SSZ_Negative_Transmit_Power_MW[ssz, tmp]
                                     * model.ssz_negative_direction_hurdle_rate_per_mw[ssz, model.period[tmp]]) \
                    * timepoint_weight

        # ### Penalties ### #
        # Unserved energy and overgeneration penalties for each zone
        for z in model.ZONES:
            # Unserved energy costs
            if model.allow_unserved_energy:
                unserved_energy_costs += model.Unserved_Energy_MW[z, tmp] * model.unserved_energy_penalty_per_mw \
                    * timepoint_weight
            # Overgen costs
            overgen_costs += model.Overgeneration_MW[z, tmp] * model.overgen_penalty_per_mw \
                * timepoint_weight

        # Reserve violation penalties
        # spinning reserve violation costs
        spin_violation_costs += model.Spin_Violation_MW[tmp] \
            * model.spin_violation_penalty_per_mw * timepoint_weight

        # Upward regulation violation costs
        upward_reg_violation_costs += model.Upward_Reg_Violation_MW[tmp] \
            * model.upward_reg_violation_penalty_per_mw * timepoint_weight
        # Downward regulation violation costs
        downward_reg_violation_costs += model.Downward_Reg_Violation_MW[tmp] \
            * model.downward_reg_violation_penalty_per_mw * timepoint_weight

        # Upward load-following reserve violation costs
        upward_lf_reserve_violation_costs += model.Upward_LF_Reserve_Violation_MW[tmp] \
            * model.upward_lf_reserve_violation_penalty_per_mw * timepoint_weight
        # Downward load-following reserve violation costs
        downward_lf_reserve_violation_costs += model.Downward_LF_Reserve_Violation_MW[tmp] \
            * model.downward_lf_reserve_violation_penalty_per_mw * timepoint_weight

    total_cost = \
        capital_and_fixed_resource_costs + \
        deliverability_cost + \
        variable_costs + \
        fuel_costs + \
        unit_start_costs + \
        unit_shutdown_costs + \
        hurdle_rate_costs + \
        unserved_energy_costs + \
        overgen_costs + \
        scheduled_curtailment_costs + \
        upward_reg_violation_costs + \
        downward_reg_violation_costs + \
        upward_lf_reserve_violation_costs + \
        downward_lf_reserve_violation_costs + \
        spin_violation_costs

    if model.allow_tx_build:
        total_cost += transmission_build_costs

    return total_cost


###########################
# ##### CONSTRAINTS ##### #
###########################

# ##### Reliability Constraints ##### #

def transmission_reliability_contribution_rule(model, reliability, period):
    """
    reliability contribution of Tx lines can come from A) Existing Tx contribution
    and B) New Tx contribution

    Args:
      model:
      period:

    Returns:

    """
    transmission_reliability_contribution = 0
    for tx_reliability, tx_line, year in model.TRANSMISSION_RELIABILITY_MAP:
        if year == period and tx_reliability == reliability:
            transmission_reliability_contribution += \
                model.Existing_Tx_Reliability_Contribution_MW[reliability, tx_line, period]
            if model.allow_tx_build:
                transmission_reliability_contribution += \
                    model.New_Tx_Reliability_Contribution_MW[reliability, tx_line, period]

    return transmission_reliability_contribution

def resource_non_elcc_surface_reliability_contribution_rule(model, reliability, period):
    non_elcc_reliability_capacity = 0.0
    for resource in model.RELIABILITY_RESOURCES[reliability]:
        if model.allow_elcc_surfaces:
            if resource not in model.ELCC_RESOURCES[reliability]:
                non_elcc_reliability_capacity += model.Reliability_Capacity_NQC_MW[reliability, resource, period]
        else:
            non_elcc_reliability_capacity += model.Reliability_Capacity_NQC_MW[reliability, resource, period]

    return non_elcc_reliability_capacity


def reliability_capacity_period_rule(model, reliability, period):
    """Ensure that there is enough capacity to meet a percentage - the reliability_margin - above peak load.

    Args:
      model:
      period:

    Returns:

    """
    # Sum the capacity of non ELCC resources,
    # This will generally include thermal, (large-scale) hydro, conventional DR, and baseload RPS resources
    # net_qualifying_capacity_fraction derates the capacity that is used in the reliability constraint
    # relative to the installed capacity of the resource
    non_elcc_surface_capacity = model.Resource_Non_ELCC_Surface_Reliability_Contribution_MW[reliability, period]

    # Sum elcc mw for all surfaces
    elcc_mw = 0.0
    if model.allow_elcc_surfaces:
        for val in model.ELCC_SURFACE_NAMES:
            if val[0]==reliability:
                elcc_mw += model.ELCC_MW[period, val[0], val[1]]

    tx_reliability_contribution = 0
    if model.allow_tx_reliability_contribution:
        tx_reliability_contribution = model.Transmission_Reliability_Contribution_MW[reliability, period]

    return non_elcc_surface_capacity \
        + tx_reliability_contribution \
        + elcc_mw \
        + model.reliability_capacity_adjustment_mw[reliability, period] \
        >= model.reliability_peak_load_mw[reliability, period] * (1 + model.reliability_margin[reliability, period])

def reliability_preNQC_power_rule(model, reliability, resource, period):
    # capacity used for NQC & ELCC derating should be less than or equal to
    if resource in model.RELIABILITY_ENERGY_LIMITED_RESOURCES:
        return (model.Reliability_Capacity_preNQC_MW[reliability, resource, period]
            <= model.Reliability_Capacity_Delivered_MW[resource, period])
    else:
        return (model.Reliability_Capacity_preNQC_MW[reliability, resource, period]
                == model.Reliability_Capacity_Delivered_MW[resource, period])


def reliability_duration_MWh_rule(model, resource, period):
    """
    first get delivered MW (derate for non-deliverable resources) for reliability MWh calc
    then calculate reliability duration MWh based on resource type

    """
    if resource in model.STORAGE_RESOURCES:
        return model.Total_Storage_Energy_Capacity_MWh[resource, period]
    elif resource in model.CONVENTIONAL_DR_RESOURCES:
        return model.conventional_dr_reliability_hours_per_day[resource, period] \
               * model.Reliability_Capacity_Delivered_MW[resource, period]


def reliability_preNQC_energy_rule(model, reliability, resource, period):
    # derates for short-duration storage & short-length shed dr resources
    """
        if resource is dr & horizon hours >24:
        preNQC_MW <= daily available capacity
        else:
        For each resource, preNQC_MW < = delivered_MWh/event window
        """
    if resource in model.RELIABILITY_ENERGY_LIMITED_RESOURCES:
        if resource in model.CONVENTIONAL_DR_RESOURCES:
            return (model.Reliability_Capacity_preNQC_MW[reliability, resource, period]
                    <= model.Reliability_Duration_MWh[resource, period] / min(24, model.reliability_event_length_hour[reliability, period]))
        else:
            return (model.Reliability_Capacity_preNQC_MW[reliability, resource, period]
                    <= model.Reliability_Duration_MWh[resource, period] / model.reliability_event_length_hour[reliability, period])
    else:
        return pyo.Constraint.Skip

def reliability_preNQC_conventional_dr_annual_rule(model, reliability, resource, period):
    # only applies to DR resources under the reliability constraint
    if resource in model.CONVENTIONAL_DR_RESOURCES:
        return (model.Reliability_Capacity_preNQC_MW[reliability, resource, period]
                <= model.Reliability_Capacity_Delivered_MW[resource, period]
                * model.conventional_dr_availability_hours_per_year[resource, period]
                / model.reliability_event_length_hour[reliability, period])
    else:
        return pyo.Constraint.Skip


def _nqc_validation_rule(model, value, reliability, resource, period):
    if resource in model.LOAD_ONLY_RESOURCES:
        return value == 0
    else:
        return True


def reliability_capacity_nqc_derate_rule(model, reliability, resource, period):
    return model.Reliability_Capacity_preNQC_MW[reliability, resource, period] \
           * model.net_qualifying_capacity_fraction[reliability, resource, period]


def elcc_surface_axes_count_rule(model, reliability_constraint, elcc_surface_name):

    axis_list = []
    for val in model.ELCC_SURFACE_MAP:
        if val[0]==reliability_constraint and val[1]==elcc_surface_name:
            axis_list.append(model.elcc_surface_axes[val])

    return max(axis_list)


def elcc_resources_init(model, reliability_constraint):
    resource_list = []
    for reliability, surface, resource in model.ELCC_SURFACE_MAP:
        if reliability==reliability_constraint:
            resource_list.append(resource)
    # make sure one resource doesn't show up on multiple surfaces within the same reliability constraint
    if len(resource_list) == len(set(resource_list)):
        return sorted(set(resource_list))
    else:
        raise ValueError(
            'Reliability {} resource assignment to ELCC surfaces is non-unique. '
            'Look for resources appearing on multiple ELCC surfaces under the reliability constraint.'
                .format(reliability_constraint))


def elcc_surface_rule(model, period, reliability_constraint, elcc_surface_name, facet):
    facet_elcc_mw = model.Elcc_Surface_Facet_LHS[period, reliability_constraint, elcc_surface_name, facet]

    return model.ELCC_MW[period, reliability_constraint, elcc_surface_name]<=facet_elcc_mw


def elcc_surface_facet_LHS_rule(model, period, reliability_constraint, elcc_surface_name, facet):
    facet_elcc_mw = 0.0
    # for each axis, add up all resources (NQC MW * multiplier) on the axis;
    # then apply axis coefficient for the facet
    # them multiply by peak MW / annual load MWh
    for axis in range(1, model.elcc_surface_axes_count[reliability_constraint, elcc_surface_name] + 1):
        axis_elcc_mw = 0.0
        for val in model.ELCC_SURFACE_MAP:
            if val[0] == reliability_constraint and val[1] == elcc_surface_name \
                    and model.elcc_surface_axes[val] == axis:
                axis_elcc_mw += model.Reliability_Capacity_NQC_MW[val[0], val[2], period] \
                                * model.elcc_surface_axes_map_multiplier[val]
        facet_elcc_mw += axis_elcc_mw \
                         * model.elcc_surface_facets_axes_coefficient[
                             reliability_constraint, elcc_surface_name, facet, axis, period]
        # add intercept
    facet_elcc_mw += model.elcc_surface_facets_axes_coefficient[
        reliability_constraint, elcc_surface_name, facet, 0, period]
    return facet_elcc_mw




# Retirements ########
# Retirement for planned and new capacity are treated separately because some constraints

# Minimum cumulative planned retirements for every period (including forced and economic planned capacity retirements)
def minimum_planned_retirements_rule(model, resource, period):
    """

    Args:
      model:
      resource:
      period:

    Returns:

    """
    # NOTE: For the first period there is no forced retirement
    # because the planned capacity in the first period (an input) should already be net of any planned retirements.
    # The model can still choose economically retire planned capacity in the first period via other constraints.
    if period == model.first_period:
        return pyo.Constraint.Skip
    else:
        # the planned capacity input values may decrease from period to period, implying retirements by a certain date
        # this constraint forces retirements by that date, but also allows for the capacity to be retired earlier
        return model.cumulative_planned_capacity_subtractions_mw[resource, period] <= \
               model.Retire_Planned_Capacity_Cumulative_MW[resource, period]


# Minimum_Planned_Retirements_Constraint
# The planned capacity in the first period is equal to the cumulative planned additions
# this constraint works in concert with Minimum_Planned_Retirements_Constraint
# to force scheduled retirements by a specific date while also allowing the model freedom to retire resources earlier
def maximum_planned_retirements_rule(model, resource, period):
    """

    Args:
      model:
      resource:
      period:

    Returns:

    """
    return model.Retire_Planned_Capacity_Cumulative_MW[resource, period] <= \
           model.cumulative_planned_capacity_additions_mw[resource, period]


def retire_new_increasing_only_rule(model, resource, period, vintage):
    """

    Args:
      model:
      resource:
      period:
      vintage:

    Returns:

    """
    # retired capacity of each vintage is strictly increasing overtime - no mothballing

    # Don't allow retirements the year in which the resource was installed
    # there should be natural disincentives to do this, but just in case explicitly remove the option.
    # Also, the else below looks back to the previous period, which wouldn't exist when period == vintage
    if period == vintage:
        return model.Retire_New_Capacity_By_Vintage_Cumulative_MW[resource, period, vintage] == 0
    # From the second period onwards, retire new capacity that was built in the previous period.
    else:
        return model.Retire_New_Capacity_By_Vintage_Cumulative_MW[resource, period, vintage] >= \
               model.Retire_New_Capacity_By_Vintage_Cumulative_MW[resource, find_prev_period(model, period), vintage]


def retire_new_upper_bound_rule(model, resource, period, vintage):
    """

    Args:
      model:
      resource:
      period:
      vintage:

    Returns:

    """
    # can't retire more capacity of a given vintage than was installed in the first place
    return model.Retire_New_Capacity_By_Vintage_Cumulative_MW[resource, period, vintage] <= \
           model.Build_Capacity_MW[resource, vintage]


def min_operational_planned_capacity_rule(model, resource, period):
    """Retain (don't retire) a minimum amount of planned capacity in each period.
    Typically used to model factors outside the scope of the model that cause generation capacity to be retained.

    Args:
      model:
      resource:
      period:

    Returns:

    """
    if model.min_operational_planned_capacity_mw[resource, period] == 0:
        return pyo.Constraint.Skip
    else:
        return model.Operational_Planned_Capacity_MW[resource, period] \
               >= model.min_operational_planned_capacity_mw[resource, period]


def min_cumulative_new_build_rule(model, new_build_resource, period):
    """For resources that can be built, how much cumulative capacity must be built by a certain year.
    Note that there is not currently an analogous constraint for storage energy capacity min new build,
    but any storage resource for which power capacity is installed via this constraint
    will be forced to have at least min_duration_h over the total power capacity of the resource.

    Args:
      model:
      period:
      new_build_resource:

    Returns:

    """
    return model.Operational_New_Capacity_MW[new_build_resource, period] \
           >= model.min_cumulative_new_build_mw[new_build_resource, period]


def capacity_limit_rule(model, capacity_limited_resource, period):
    """Limit on amount of capacity that can be installed for capacity-limited resources.
    For energy efficiency the limit should be specified in average MWs (aMW)

    Args:
      model:
      period:
      capacity_limited_resource:

    Returns:

    """
    return model.Operational_Capacity_MW[capacity_limited_resource, period] \
           <= model.capacity_limit_mw[capacity_limited_resource, period]


def custom_constraint_period_lhs(model, constraint_group, period):
    """
    The expression makes the LHS of the generic custom_constraint_period_rule.
    index_2 allows for variables and expressions based on memberships between objects to be included
    Many variables and params are defined on only one object/index - in this case set index_2=None in the input data.
    Args:
      model:
      constraint_group:
      period:

    Returns:

    """
    lhs_sum = 0.0

    for group, index_1, index_2, attribute in model.CUSTOM_CONSTRAINT_PERIOD_GROUP_MAP:
        if group == constraint_group:
            if index_2 == 'None':
                lhs_sum += getattr(model, attribute)[index_1, period] * \
                           model.custom_constraint_period_multiplier[group, index_1, index_2, attribute]
            else:
                lhs_sum += getattr(model, attribute)[index_1, index_2, period] * \
                           model.custom_constraint_period_multiplier[group, index_1, index_2, attribute]

    return lhs_sum


def custom_constraint_period_rule(model, constraint_group, period):
    """

    Args:
      model:
      constraint_group:
      period:

    Returns:

    """
    if model.custom_constraint_period_limit[constraint_group, period] == float('inf'):
        return pyo.Constraint.Skip
    else:
        return (model.Custom_Constraint_Period_LHS[constraint_group, period] <=
                model.custom_constraint_period_limit[constraint_group, period]
                )


def custom_constraint_timepoint_lhs(model, constraint_group, timepoint):
    """
    The expression makes the LHS of the generic custom_constraint_timepoint_rule.
    index_2 allows for variables and expressions based on memberships between objects to be included
    Many variables and params are defined on only one object/index - in this case set index_2=None in the input data.

    Args:
      model:
      constraint_group:
      timepoint:

    Returns:

    """
    lhs_sum = 0.0

    for group, index_1, index_2, attribute in model.CUSTOM_CONSTRAINT_TIMEPOINT_GROUP_MAP:
        if group == constraint_group:
            if index_2 == 'None':
                lhs_sum += getattr(model, attribute)[index_1, timepoint] * \
                           model.custom_constraint_timepoint_multiplier[group, index_1, index_2, attribute]
            else:
                lhs_sum += getattr(model, attribute)[index_1, index_2, timepoint] * \
                           model.custom_constraint_timepoint_multiplier[group, index_1, index_2, attribute]

    return lhs_sum


def custom_constraint_timepoint_rule(model, constraint_group, timepoint):
    """

    Args:
      model:
      constraint_group:
      timepoint:

    Returns:

    """
    if model.custom_constraint_timepoint_limit[constraint_group, timepoint] == float('inf'):
        return pyo.Constraint.Skip
    else:
        return (model.Custom_Constraint_Timepoint_LHS[constraint_group, timepoint] <=
                model.custom_constraint_timepoint_limit[constraint_group, timepoint]
                )


def storage_minimum_duration_rule(model, new_storage_resource, period):
    """Minimum duration for new storage resources.

    Args:
      model:
      period:
      new_storage_resource:

    Returns:

    """
    return model.Cumulative_New_Storage_Energy_Capacity_MWh[new_storage_resource, period] \
           >= model.Operational_New_Capacity_MW[new_storage_resource, period] \
           * model.min_duration_h[new_storage_resource]


# ## Fuel ## #
def fuel_consumption_rule(model, resource, timepoint):
    """Constraint that represents the fuel consumption of each thermal resource in each timepoint.
    For resources without a blend fuel,
    this constraint is essentially an expression that defines Base_Fuel_Consumption_MMBtu.
    For resource with a blend fuel, it forces the model to choose between base and blend fuels.
    Includes the dispatch up/down from providing reserves and the associated fuel increase/savings.

    Args:
      model:
      resource:
      timepoint:

    Returns:

    """

    power_output_mwh = model.Provide_Power_MW[resource, timepoint]

    # Regulation reserve dispatch
    if resource in model.REGULATION_RESERVE_RESOURCES:
        power_output_mwh += (model.Provide_Upward_Reg_Net_MW[resource, timepoint]
                             - model.Provide_Downward_Reg_Net_MW[resource, timepoint]) \
                             * model.reg_dispatch_fraction

    # Load-following reserve dispatch
    if resource in model.LOAD_FOLLOWING_RESERVE_RESOURCES:
        power_output_mwh += (model.Provide_LF_Upward_Reserve_Net_MW[resource, timepoint]
                             - model.Provide_LF_Downward_Reserve_Net_MW[resource, timepoint]) \
                             * model.lf_reserve_dispatch_fraction

    # Calculate fuel burn using the fuel burn slope
    fuel_consumption = power_output_mwh \
        * model.fuel_burn_slope_mmbtu_per_mwh[resource]

    # add intercept times commitment for unit commitment resources
    if resource in model.UNIT_COMMITMENT_RESOURCES:
        fuel_consumption += model.Commit_Units[resource, timepoint] \
            * model.fuel_burn_intercept_mmbtu_per_hr[resource]

        # add start fuel burn
        fuel_consumption += model.Start_Fuel_MMBtu[resource, timepoint]

    # if the resource' operation allow blending fuel
    if resource in model.BLEND_FUEL_RESOURCES:
        return fuel_consumption == (model.Base_Fuel_Consumption_MMBtu[resource, timepoint]
                                    + model.Blend_Fuel_Consumption_MMBtu[resource, timepoint])
    else:
        return fuel_consumption == model.Base_Fuel_Consumption_MMBtu[resource, timepoint]


# ##### RPS constraints ##### #
def define_previously_banked_rps(model, policy, period):
    """RPS credits banked from previous period or planned spending of historical bank

    Args:
      model:
      policy:
      period:

    Returns:

    """
    if not model.optimize_rps_banking:
        # if optimized banking is not enabled, no rps credits are banked from any previous period
        # Legacy: (Replacing planned_spend) If users wish to specify planned spending of RPS bank,
        # turn off optimize rps_banking, and use the non_modeled_RPS input by period
        return 0
    else:
        if period == model.first_period:
            # the first period starts with starting_rps_bank_mwh
            # starting_rps_bank_mwh can be zero if you don't want to start with any bank
            return model.starting_rps_bank_mwh[policy] / model.years_in_period[period]
        else:
            # Each period can have a different length (years_in_period) so rps credits that are banked
            # in a period must be multiplied by the years_in_period when the credits go into the bank (in prev_period).
            # The credits are divided by years_in_period when they come out of the bank in the subsequent period
            prev_period = find_prev_period(model, period)
            return model.Bank_RPS_MWh[policy, prev_period] * \
                model.years_in_period[prev_period] / model.years_in_period[period]


def define_yearly_rps_storage_losses(model, policy, period):
    """Derived variable that represents the total amount of storage losses in RPS zones over the course of a year.
    This will be zero if storage losses aren't counted against the RPS requirement.

    Storage losses count as curtailment on the assumption that RPS-eligible energy is going into storage.
    Storage losses are incurred when providing reserves because reserves will be dispatched some fraction of the time:
    the reserve_dispatch_fraction.

    Args:
      model:
      policy:
      period:

    Returns:

    """
    yearly_storage_losses = 0.0

    # skip this constraint when storage loss
    if not model.count_storage_losses_as_rps_curtailment:
        return pyo.Constraint.Skip

    # if resources are
    for storage_resource in model.STORAGE_RESOURCES:
        if model.zone[storage_resource] in model.RPS_ZONES[policy]:
            for timepoint in model.TIMEPOINTS:
                if model.period[timepoint] == period:

                    charging_mwh = 0.0
                    discharging_mwh = 0.0

                    charging_mwh += model.Increase_Load_MW[storage_resource, timepoint]
                    discharging_mwh += model.Provide_Power_MW[storage_resource, timepoint]

                    # Regulation reserve dispatch
                    if storage_resource in model.REGULATION_RESERVE_RESOURCES:
                        charging_mwh += model.Provide_Downward_Reg_From_Charge_MW[storage_resource, timepoint] \
                            * model.reg_dispatch_fraction
                        discharging_mwh += model.Provide_Upward_Reg_MW[storage_resource, timepoint] \
                            * model.reg_dispatch_fraction

                    # Load-following reserve dispatch
                    if storage_resource in model.LOAD_FOLLOWING_RESERVE_RESOURCES:
                        charging_mwh += model.Provide_LF_Downward_Reserve_From_Charge_MW[storage_resource, timepoint] \
                            * model.lf_reserve_dispatch_fraction
                        discharging_mwh += model.Provide_LF_Upward_Reserve_MW[storage_resource, timepoint] \
                            * model.lf_reserve_dispatch_fraction

                    # Calculate losses
                    # To calculate a yearly total, losses are weighted by day_weight
                    # The values charging_mwh/discharging_mwh represent the load/generation as seen by the grid
                    # For example, discharging_mwh represents energy put onto the grid by the storage device,
                    # which is less than the amount of energy taken
                    # out of the storage device if discharging_efficiency is < 100%
                    yearly_storage_losses += \
                        charging_mwh * model.day_weight[model.day[timepoint]] \
                        * (1.0 - model.charging_efficiency[storage_resource]) \
                        + discharging_mwh * model.day_weight[model.day[timepoint]] \
                        * (1.0 / model.discharging_efficiency[storage_resource] - 1.0)

    return yearly_storage_losses


def define_yearly_rps_eligible_generation(model, policy, period):
    """How much renewable power contracted to all zones included in the RPS is produced in each period?
    This derived variable represents the amount of power actually produced by generators,
    not the pre-curtailment potential.

    Args:
      model:
      policy:
      period:

    Returns:

    """

    total_rps_eligible_generation = 0.0

    for rps_r in model.RPS_REGULAR_RESOURCES[policy]:
        for timepoint in model.TIMEPOINTS:
            if model.period[timepoint] == period:
                total_rps_eligible_generation += (model.Provide_Power_MW[rps_r, timepoint]
                                                 * model.day_weight[model.day[timepoint]]
                                                 * model.rps_eligible_fraction[policy, rps_r])

    return total_rps_eligible_generation


def define_scheduled_rps_curtailment_per_year(model, policy, period):
    """Derived variable that represents the total MWh curtailed from RPS-eligible resources.

    Args:
      model:
      policy:
      period:

    Returns:

    """
    scheduled_rps_curtailment = 0.0

    for resource in model.RPS_ELIGIBLE_RESOURCES[policy]:
        if resource in model.CURTAILABLE_VARIABLE_RESOURCES:
            for timepoint in model.TIMEPOINTS:
                if model.period[timepoint] == period:
                    scheduled_rps_curtailment += model.Scheduled_Curtailment_MW[resource, timepoint] \
                                                 * model.day_weight[model.day[timepoint]]* model.rps_eligible_fraction[policy, resource]

    return scheduled_rps_curtailment


def define_subhourly_rps_curtailment_per_year(model, policy, period):
    """Derived variable that represents the total MWh curtailed from RPS-eligible resources on the sub-hourly timescale.

    Args:
      model:
      policy:
      period:

    Returns:

    """
    total_subhourly_curtailment = 0.0

    for timepoint in model.TIMEPOINTS:
        if model.period[timepoint] == period:
            for zone in model.RPS_ZONES[policy]:
                if zone in model.LOAD_FOLLOWING_ZONES:
                    total_subhourly_curtailment += (model.day_weight[model.day[timepoint]]
                                                    *(model.Subhourly_Downward_LF_Energy_MWh[zone, timepoint]
                                                      - model.Subhourly_Upward_LF_Energy_MWh[zone, timepoint]))

    return total_subhourly_curtailment


def achieve_rps_rule(model, policy, period):
    """Ensure that the RPS target is met in each period.

    Args:
      model:
      period:

    Returns:

    """

    # Bank_RPS_MWh is the slack variable - if optimize_rps_banking is false, then Bank_RPS_MWh will still represent
    # any additional rps credits that are generated above the rps target, but these credits are not allowed to count
    # towards the rps target in the subsequent period
    # because in this case Previously_Banked_RPS_MWh does not include the variable Bank_RPS_MWh[prev_period]
    return model.Previously_Banked_RPS_MWh[policy, period] \
        + model.RPS_Eligible_Generation_MWh_Per_Year[policy, period] \
        + model.RPS_Blend_Fuel_Generation_MWh_Per_Year[policy, period] \
        - model.RPS_Storage_Losses_MWh_Per_Year[policy, period] \
        - model.Subhourly_RPS_Curtailment_Per_Year[policy, period] \
        + model.rps_nonmodeled_mwh[policy, period] \
        == model.RPS_Target_MWh[policy, period] + model.Bank_RPS_MWh[policy, period]


def define_unbundled_flag(model, policy, period):
    """Determine whether the model will enforce the RPS_Unbundled_Fraction_Limit_Constraint

    Args:
      model:
      policy:
      period:

    Returns:

    """
    number_of_unbundled_resources = 0
    for rps_r in model.RPS_ELIGIBLE_RESOURCES[policy]:
        if model.zone[rps_r] not in model.RPS_ZONES[policy]:
            number_of_unbundled_resources += 1

    # Don't write constraint if there aren't unbundled resources.
    if number_of_unbundled_resources == 0:
        return False
    # Don't write constraint if an RPS target isn't enforced in the period.
    elif model.rps_fraction_of_retail_sales[policy, period] == 0:
        return False
    else:
        return True


def define_yearly_unbundled_rps_generation(model, policy, period):
    """An expression that calculates the amount of unbundled rps for each period
    Unbundled here is equivalent to out of zones, i.e. How much renewable power contracted to RPS zones but balanced
    elsewhere is produced in each period?

    :param model:
    :param policy:
    :param period:
    :return:
    """

    yearly_unbundled_rps_generation = 0.0

    for rps_r in model.RPS_ELIGIBLE_RESOURCES[policy]:
        if (model.zone[rps_r] not in model.RPS_ZONES[policy]) and (rps_r not in model.EE_PROGRAMS):
            for timepoint in model.TIMEPOINTS:
                if model.period[timepoint] == period:
                    # Sum the amount of power actually produced by generators (not the pre-curtailment potential).
                    yearly_unbundled_rps_generation += (model.Provide_Power_MW[rps_r, timepoint] *
                                                        model.day_weight[model.day[timepoint]]*
                                                        model.rps_eligible_fraction[policy, rps_r])

    return yearly_unbundled_rps_generation


def rps_unbundled_fraction_limit_rule(model, policy, period):
    """Limit the amount of unbundled RECs to a fraction of the total RPS target.
    If rps_unbundled_fraction_limit = 1 then the limit isn't binding.
    If rps_unbundled_fraction_limit = 0 then no unbundled RECs are allowed.
    Be careful with = 0 as this could cause an infeasibility with existing/planned unbundled non-curtailable resources
    or curtailment of all power from existing/planned unbundled curtailable resources.
    Storage losses and subhourly curtailment in the location zone of the unbundled resources
    aren't currently included in this limit.

    Args:
      model:
      policy:
      period:

    Returns:

    """

    if model.enforce_unbundled_fraction_limit[policy, period]:
        return model.Unbundled_RPS_Generation_MWh_Per_Year[policy, period] <=\
               (model.rps_unbundled_fraction_limit[policy, period] * model.RPS_Target_MWh[policy, period])
    else:
        return pyo.Constraint.Skip


# ##### Emissions Target Constraints ##### #
# Emissions constraint by period
def emissions_target_rule(model, group, period):
    """If an emissions limit is enforced for the emissions group, the emissions from that group must be less than or equal to a period-varying target.
    A "credit" is included that can be either positive or negative to model adjustments to the target.

    Args:
      model:
      group_constraint_period:

    Returns:

    """
    return model.Resource_Emissions_Units_Per_Year[group, period] \
           + model.Unspecified_Imports_Units_Per_Year[group, period] <= \
           model.emissions_target_units_per_year[group, period] \
           + model.emissions_credit_units_per_year[group, period]


# ##### Resource Operational Constraints ##### #


# ### Variable renewables ### #
def variable_power_rule(model, variable_resource, timepoint):
    """The amount of power that is produced by variable renewable resources is equal to the potential production (shape)
    multiplied by an available_capacity_fraction if applicable,
    If a resource can be scheduled to be curtailed, then a curtailment variable is included.
    Note: subhourly curtailment is treated elsewhere - this constraint only includes scheduled curtailment.
    The model could have numeric problems with this constraint when model.shape is very small,
    so if the production potential is very small, round to zero via shape_trimmed

    Args:
      model:
      timepoint:
      variable_resource:

    Returns:
      Renewable_Power_Constraint rule

    """
    shape = model.shape[variable_resource, model.day[timepoint], model.hour_of_day[timepoint]]
    if shape < 0.0001:
        shape = 0

    if variable_resource in model.CURTAILABLE_VARIABLE_RESOURCES:
        curtail_or_not = model.Scheduled_Curtailment_MW[variable_resource, timepoint]
    else:
        curtail_or_not = 0

    return model.Provide_Power_MW[variable_resource, timepoint] + curtail_or_not == \
           model.Available_Capacity_In_Timepoint_MW[variable_resource, timepoint] \
           * shape


# ### Thermal ### #


def unit_commitment_commitment_rule(model, resource, timepoint):
    """How much capacity of each unit commitment resource can be committed (or started) in each timepoint.

    If the Param `commit_all_capacity` is defined for a resource, it will be forced on.
    Otherwise, the capacity that is committed or available to be started is bounded by the total available capacity.

    Args:
        model ([type]): [description]
        resource ([type]): [description]
        timepoint ([type]): [description]

    Returns:
        [type]: [description]
    """

    if model.commit_all_capacity[resource, timepoint]:
        return (
                model.Commit_Capacity_MW[resource,
                                         timepoint] == model.Available_Capacity_In_Timepoint_MW[resource,
                                                                                                timepoint]
        )
    else:
        return (
                model.Commit_Capacity_MW[resource,
                                         timepoint] +
                model.Starting_Units[resource,
                                     timepoint] * model.unit_size_mw[resource] <=
                model.Available_Capacity_In_Timepoint_MW[resource,
                                                         timepoint]
        )


def unit_commitment_max_gen_and_up_reserves_rule(model, resource, timepoint):
    """Maximum output and, for resources that can provide them, upward reserves and frequency response
    from a unit commitment resource:
    the sum of generation and reserves cannot exceed the fully operational (or committed) capacity
    for each resource in each timepoint.

    Args:
      model:
      timepoint:
      resource:

    Returns:
      Unit_Commitment_Max_Gen_Up_Reserve_Constraint

    """
    if resource in model.REGULATION_RESERVE_RESOURCES:
        upward_reg = model.Provide_Upward_Reg_Net_MW[resource, timepoint]
    else:
        upward_reg = 0

    if resource in model.LOAD_FOLLOWING_RESERVE_RESOURCES:
        upward_lf_reserves = model.Provide_LF_Upward_Reserve_Net_MW[resource, timepoint]
    else:
        upward_lf_reserves = 0

    if resource in model.SPINNING_RESERVE_RESOURCES:
        spin = model.Provide_Spin_Net_MW[resource, timepoint]
    else:
        spin = 0

    if resource in model.TOTAL_FREQ_RESP_RESOURCES:
        frequency_response = model.Provide_Frequency_Response_Net_MW[resource, timepoint]
    else:
        frequency_response = 0

    if resource in model.UNIT_COMMITMENT_RAMP_LIMITED_RESOURCES:
        max_power = \
            (model.Commit_Units[resource, timepoint] - model.Fully_Operational_Units[resource, timepoint]) \
            * model.unit_size_mw[resource] \
            * model.min_stable_level_fraction[resource] \
            + model.Fully_Operational_Units[resource, timepoint] \
            * model.unit_size_mw[resource]
    else:
        max_power = model.Commit_Capacity_MW[resource, timepoint]

    return model.Provide_Power_MW[resource, timepoint] \
           + upward_reg \
           + upward_lf_reserves \
           + spin \
           + frequency_response \
           <= max_power


def unit_commitment_upward_reserve_ramp_rule(model, resource, timepoint):
    """Maximum upward reserves from a unit commitment resource within the timeframe assumed
    for reg, load following, and spin reserves (reserve_timeframe_fraction_of_hour):
    The sum of reserves cannot exceed the operational/committed capacity for each resource in each timepoint.
    Frequency response reserve is not included here because it is on a faster timescale (~1 min).

    Args:
      model:
      timepoint:
      resource:

    Returns:
      Unit_Commitment_Upward_Reserve_Ramp_Constraint

    """
    if resource in model.REGULATION_RESERVE_RESOURCES:
        upward_reg = model.Provide_Upward_Reg_Net_MW[resource, timepoint]
    else:
        upward_reg = 0

    if resource in model.LOAD_FOLLOWING_RESERVE_RESOURCES:
        upward_lf_reserves = model.Provide_LF_Upward_Reserve_Net_MW[resource, timepoint]
    else:
        upward_lf_reserves = 0

    if resource in model.SPINNING_RESERVE_RESOURCES:
        spin = model.Provide_Spin_Net_MW[resource, timepoint]
    else:
        spin = 0

    if resource in model.UNIT_COMMITMENT_RAMP_LIMITED_RESOURCES:
        reserve_capable_units = model.Fully_Operational_Units[resource, timepoint]
    else:
        reserve_capable_units = model.Commit_Units[resource, timepoint]

    # don't compile constraint if unit can ramp the entire operational range (1-PMin) within the reserve timeframe
    # or if the unit can't provide any reserves
    if (
            resource not in model.REGULATION_RESERVE_RESOURCES and
            resource not in model.LOAD_FOLLOWING_RESERVE_RESOURCES and resource not in model.SPINNING_RESERVE_RESOURCES
    ):
        return pyo.Constraint.Skip
    elif model.reserve_timeframe_fraction_of_hour * model.ramp_rate_fraction[resource] >= \
            (1 - model.min_stable_level_fraction[resource]):
        return pyo.Constraint.Skip
    else:
        return upward_reg + upward_lf_reserves + spin <= \
               reserve_capable_units \
               * model.unit_size_mw[resource] \
               * model.ramp_rate_fraction[resource] \
               * model.reserve_timeframe_fraction_of_hour


def unit_commitment_downward_reserve_ramp_rule(model, resource, timepoint):
    """Maximum downward reserves from a unit commitment resource within the timeframe assumed
    for reg and load following reserves (reserve_timeframe_fraction_of_hour):
    The sum of reserves cannot exceed the operational/committed capacity for each resource in each timepoint.

    Args:
      model:
      timepoint:
      resource:

    Returns:
      Unit_Commitment_Downward_Reserve_Ramp_Constraint

    """
    if resource in model.REGULATION_RESERVE_RESOURCES:
        downward_reg = model.Provide_Downward_Reg_Net_MW[resource, timepoint]
    else:
        downward_reg = 0

    if resource in model.LOAD_FOLLOWING_RESERVE_RESOURCES:
        downward_lf_reserves = model.Provide_LF_Downward_Reserve_Net_MW[resource, timepoint]
    else:
        downward_lf_reserves = 0

    if resource in model.UNIT_COMMITMENT_RAMP_LIMITED_RESOURCES:
        reserve_capable_units = model.Fully_Operational_Units[resource, timepoint]
    else:
        reserve_capable_units = model.Commit_Units[resource, timepoint]

    # don't compile constraint if unit can ramp the entire operational range (1-PMin) within the reserve timeframe
    # or if the unit can't provide any reserves
    if resource not in model.REGULATION_RESERVE_RESOURCES and resource not in model.LOAD_FOLLOWING_RESERVE_RESOURCES:
        return pyo.Constraint.Skip
    elif model.reserve_timeframe_fraction_of_hour * model.ramp_rate_fraction[resource] >= \
            (1 - model.min_stable_level_fraction[resource]):
        return pyo.Constraint.Skip
    else:
        return downward_reg + downward_lf_reserves <= \
               reserve_capable_units \
               * model.unit_size_mw[resource] \
               * model.ramp_rate_fraction[resource] \
               * model.reserve_timeframe_fraction_of_hour


def unit_commitment_freq_resp_limit(model, resource, timepoint):
    """Each MW of online thermal units can provide a fraction of its committed capacity of frequency response.
    (the fraction is the param thermal_freq_response_fraction_of_commitment)
    Units that are in the process of starting or stopping are assumed to not provide frequency response.

    Args:
      model:
      timepoint:
      resource:

    Returns:

    """
    if resource in model.UNIT_COMMITMENT_RAMP_LIMITED_RESOURCES:
        frequency_responsive_units = model.Fully_Operational_Units[resource, timepoint]
    else:
        frequency_responsive_units = model.Commit_Units[resource, timepoint]

    return model.Provide_Frequency_Response_Net_MW[resource, timepoint] \
           <= model.thermal_freq_response_fraction_of_commitment[resource] \
           * model.unit_size_mw[resource] \
           * frequency_responsive_units


def unit_commitment_min_gen_and_down_reserves_rule(model, resource, timepoint):
    """Minimum output and downward reserves from a unit commitment resource.

    Args:
      model:
      timepoint:
      resource:

    Returns:
      Min_Gen_Downward_Reserve_Constraint

    """
    if resource in model.REGULATION_RESERVE_RESOURCES:
        downward_reg = model.Provide_Downward_Reg_Net_MW[resource, timepoint]
    else:
        downward_reg = 0

    if resource in model.LOAD_FOLLOWING_RESERVE_RESOURCES:
        downward_lf_reserves = model.Provide_LF_Downward_Reserve_Net_MW[resource, timepoint]
    else:
        downward_lf_reserves = 0

    # Relax the min gen constraint by a tiny amount if the unit is ramp limited and is starting or stopping
    # in order to avoid numerical conflicts between the ramp up, ramp down, and maximum power constraints
    if resource in model.UNIT_COMMITMENT_RAMP_LIMITED_RESOURCES:
        ramp_relax_start_stop = model.ramp_relax \
                                * (model.Commit_Units[resource, timepoint] - model.Fully_Operational_Units[
            resource, timepoint])
    else:
        ramp_relax_start_stop = 0

    return model.Provide_Power_MW[resource, timepoint] - downward_reg - downward_lf_reserves \
           >= model.Commit_Capacity_MW[resource, timepoint] \
           * model.min_stable_level_fraction[resource] \
           - ramp_relax_start_stop


def unit_commitment_start_rule(model, resource, timepoint):
    """Can't prestart more capacity than is off (not committed, not already starting).

    Args:
      model:
      timepoint:
      resource:

    Returns:

    """
    return model.PreStart_Units[resource, timepoint] * model.unit_size_mw[resource] \
           + model.Starting_Units[resource, timepoint] * model.unit_size_mw[resource] \
           + model.Commit_Capacity_MW[resource, timepoint] \
           <= model.Available_Capacity_In_Timepoint_MW[resource, timepoint]


def unit_commitment_shut_down_rule(model, unit_commitment_resource, timepoint):
    """Can't preshut down more units than are on (committed) and not already scheduled to shut down.

    Args:
      model:
      timepoint:

    Returns:

    """
    return model.PreShut_Down_Units[unit_commitment_resource, timepoint] \
        <= model.Commit_Units[unit_commitment_resource, timepoint] \
        - model.Shutting_Down_Units[unit_commitment_resource, timepoint]


def unit_commitment_resource_commitment_tracking_rule(model, resource, timepoint):
    """Track how many units are committed/on: the number of committed units in each timepoint is the number that were
    committed in the previous timepoint plus all units that are scheduled to start in the current timepoint minus
    all units that are scheduled to shut down in the current timepoint.

    Args:
      model:
      timepoint:
      resource:

    Returns:

    """
    return model.Commit_Units[resource, timepoint] \
           == model.Commit_Units[resource, model.previous_timepoint[timepoint]] \
           + model.Start_Units[resource, timepoint] \
           - model.Shut_Down_Units[resource, timepoint]


def unit_commitment_resource_ramp_up_rule(model, resource, timepoint):
    """Sets the upper bound ramp for a resource.
    When starting a ramp-limited resource the unit starts to PMin in the first timestep,
    so the full operational range is not available.
    ramp_relax is a small number that ensures that the ramp constraints don't conflict
    with the pmax or pmin constraints when turning on or off generators.
    Don't compile ramp rate constraints if the unit can ramp its entire operable range (1-PMin) in the timestep.
    Note: current formulation assumes hourly timesteps

    Args:
      model:
      timepoint:
      resource:

    Returns:

    """
    if model.ramp_rate_fraction[resource] >= \
            (1 - model.min_stable_level_fraction[resource]):
        return pyo.Constraint.Skip
    else:
        return model.Provide_Power_MW[resource, timepoint] \
               - model.Provide_Power_MW[resource, model.previous_timepoint[timepoint]] \
               <= (model.Commit_Capacity_MW[resource, timepoint] - model.Start_Capacity_MW[resource, timepoint]) \
               * model.ramp_rate_fraction[resource] \
               + model.Start_Capacity_MW[resource, timepoint] \
               * model.min_stable_level_fraction[resource] \
               + model.Start_Units[resource, timepoint] * model.ramp_relax


def unit_commitment_resource_ramp_down_rule(model, resource, timepoint):
    """Sets the lower bound ramp for a resource.
    When stopping a ramp-limited resource, the unit held at PMin in the timestep before the Shut_Down_Units is non-zero.
    To be able to shut the unit off, it needs to ramp down in the current timestep by Shut_Down_Units * size * PMin.
    ramp_relax is a small number that ensures that the ramp constraints don't conflict
    with the pmax or pmin constraints when turning on or off generators.
    Don't compile ramp rate constraints if the unit can ramp its entire operable range (1-PMin) in the timestep.
    Note: current formulation assumes hourly timesteps

    Args:
      model:
      timepoint:
      resource:

    Returns:

    """
    if model.ramp_rate_fraction[resource] >= \
            (1 - model.min_stable_level_fraction[resource]):
        return pyo.Constraint.Skip
    else:
        return -(model.Provide_Power_MW[resource, timepoint]
                 - model.Provide_Power_MW[resource, model.previous_timepoint[timepoint]]) \
               <= (model.Commit_Capacity_MW[resource, timepoint] - model.Start_Capacity_MW[resource, timepoint]) \
               * model.ramp_rate_fraction[resource] \
               + model.Shut_Down_Units[resource, timepoint] \
               * model.unit_size_mw[resource] \
               * model.min_stable_level_fraction[resource] \
               + model.Shut_Down_Units[resource, timepoint] * model.ramp_relax


def generate_at_max(model, resource, timepoint):
    """Generate at max resources always produce power at a level equal to their available capacity.

    Args:
      model:
      timepoint:
      resource:

    Returns:

    """

    return model.Provide_Power_MW[resource, timepoint] \
           == model.Available_Capacity_In_Timepoint_MW[resource, timepoint]



# ### Non-unit commitment ### #
def _fraction_or_undefined_validate(model, value, resource, index):
    """Validate that input param value is a fraction, or equal to "undefined"."""

    if value == 'undefined':
        return True
    else:
        return 0.0 <= value <= 1.0

def _fraction_or_undefined_validate_no_resource_that_increase_load(model, value, resource, index):
    """Validate that input param value is a fraction, or equal to "undefined"."""

    if resource in model.RESOURCES_THAT_INCREASE_LOAD:
        return value == 'undefined'
    else:
        if value == 'undefined':
            return True
        else:
            return 0.0 <=value <=1.0

def _adjacency_h_validation_rule(model, value, resource):
    if value == 9999:
        return True
    else:
        return value < 12

def energy_budget_rule(model, resource, period, day):
    """
    Generators must exactly meet a pre-defined energy budget, unless
    producing less than energy budget (aka spill, in the case of hydro) is allowed.
    Adjusts the daily energy budget for mileage when providing reserves/regulation.
    """
    # get list of timepoints for current day
    timepoints_on_day = list()
    for tmp in model.TIMEPOINTS:
        if model.period[tmp] == period and model.day[tmp] == day:
            timepoints_on_day.append(tmp)

    available_energy_current_day_mwh = (
            model.Operational_Capacity_MW[resource,
                                          period] * model.energy_budget_cap_factor[resource,
                                                                                   day] * model.timepoints_per_day
    )

    # Daily_Energy_Budget_Change_MWh can be positive or negative,
    # with negative representing a decrease in budget for the day
    if model.interday_energy_sharing:
        if resource in model.RESOURCES_WITH_MULTIDAY_SHARING:
            available_energy_current_day_mwh += model.Daily_Energy_Budget_Change_MWh[resource, period, day]

    energy_used_current_day_mwh = 0.0
    for tmp in timepoints_on_day:
        if resource in model.REGULATION_RESERVE_RESOURCES:
            upward_reg_mw = model.Provide_Upward_Reg_Net_MW[resource, tmp]
            downward_reg_mw = model.Provide_Downward_Reg_Net_MW[resource, tmp]
        else:
            upward_reg_mw = 0
            downward_reg_mw = 0

        if resource in model.LOAD_FOLLOWING_RESERVE_RESOURCES:
            upward_lf_reserves_mw = model.Provide_LF_Upward_Reserve_Net_MW[resource, tmp]
            downward_lf_reserves_mw = model.Provide_LF_Downward_Reserve_Net_MW[resource, tmp]
        else:
            upward_lf_reserves_mw = 0
            downward_lf_reserves_mw = 0

        energy_used_current_day_mwh += model.Provide_Power_MW[resource, tmp] \
                                       + (upward_reg_mw - downward_reg_mw) \
                                       * model.reg_dispatch_fraction \
                                       + (upward_lf_reserves_mw - downward_lf_reserves_mw) \
                                       * model.lf_reserve_dispatch_fraction

    if model.can_produce_less_than_energy_budget[resource]:
        return energy_used_current_day_mwh <= available_energy_current_day_mwh
    else:
        return energy_used_current_day_mwh == available_energy_current_day_mwh



def non_unit_commitment_max_gen_and_up_reserves_rule(model, resource, timepoint):
    """
    Output + up reserves cannot exceed pre-specified daily max

    For storage resources, this constraint bounds up reserves from discharge only.
    """
    if resource in model.REGULATION_RESERVE_RESOURCES:
        upward_reg = model.Provide_Upward_Reg_MW[resource, timepoint]
    else:
        upward_reg = 0

    if resource in model.LOAD_FOLLOWING_RESERVE_RESOURCES:
        upward_lf_reserves = model.Provide_LF_Upward_Reserve_MW[resource, timepoint]
    else:
        upward_lf_reserves = 0

    if resource in model.SPINNING_RESERVE_RESOURCES:
        spin = model.Provide_Spin_MW[resource, timepoint]
    else:
        spin = 0

    if resource in model.TOTAL_FREQ_RESP_RESOURCES:
        frequency_response = model.Provide_Frequency_Response_MW[resource, timepoint]
    else:
        frequency_response = 0

    return model.Provide_Power_MW[resource, timepoint] + upward_reg + upward_lf_reserves + spin + frequency_response \
           <= model.Available_Capacity_In_Timepoint_MW[resource, timepoint] * model.max_provide_power_fraction[resource, timepoint]


def min_increase_load_and_up_reserves_rule(model, resource, timepoint):
    """
    This rule forces Increase_Load_MW to be at or above Min_Increase_Load_In_Timepoint_MW.
    Any upward reserve commitment from charging counts against Increase_Load_MW because proving upward reserves
    from charging equates to potential reduction in load.
    Thus the implementation here is conservative because by including reserves in the constraint,
    the load will always be at least Min Load, even if the reserves are called upon.
    Frequently Min_Increase_Load_In_Timepoint_MW will be zero - in that situation,
    up reserves from charging are constrained to be <= Increase_Load_MW."
    """

    if resource in model.REGULATION_RESERVE_RESOURCES:
        upward_reg = model.Provide_Upward_Reg_From_Charge_MW[resource, timepoint]
    else:
        upward_reg = 0

    if resource in model.LOAD_FOLLOWING_RESERVE_RESOURCES:
        upward_lf_reserves = model.Provide_LF_Upward_Reserve_From_Charge_MW[resource, timepoint]
    else:
        upward_lf_reserves = 0

    if resource in model.SPINNING_RESERVE_RESOURCES:
        spin = model.Provide_Spin_From_Charge_MW[resource, timepoint]
    else:
        spin = 0

    if resource in model.TOTAL_FREQ_RESP_RESOURCES:
        frequency_response = model.Provide_Frequency_Response_From_Charge_MW[resource, timepoint]
    else:
        frequency_response = 0

    return model.Increase_Load_MW[resource, timepoint] - upward_reg - upward_lf_reserves - spin - frequency_response \
           >=  model.Min_Increase_Load_In_Timepoint_MW[resource, timepoint]



def non_unit_commitment_min_gen_and_down_reserves_rule(model, resource, timepoint):
    """
    Power output minus downward reserves cannot be below the pre-specified daily min gen

    For storage resources, this constraint bounds down reserves from discharge only.
    """
    if resource in model.REGULATION_RESERVE_RESOURCES:
        downward_reg = model.Provide_Downward_Reg_MW[resource, timepoint]
    else:
        downward_reg = 0

    if resource in model.LOAD_FOLLOWING_RESERVE_RESOURCES:
        downward_lf_reserves = model.Provide_LF_Downward_Reserve_MW[resource, timepoint]
    else:
        downward_lf_reserves = 0

    return model.Provide_Power_MW[resource, timepoint] - downward_reg - downward_lf_reserves \
           >= model.Available_Capacity_In_Timepoint_MW[resource, timepoint] \
           * model.min_gen_fraction[resource, timepoint]



def max_increase_load_and_down_reserves_rule(model, resource, timepoint):
    """
        Down reserves from charge + charge <= max increase load (replaces the old storage down reserve power rule)
        """
    if resource in model.REGULATION_RESERVE_RESOURCES:
        downward_reg = model.Provide_Downward_Reg_From_Charge_MW[resource, timepoint]
    else:
        downward_reg = 0

    if resource in model.LOAD_FOLLOWING_RESERVE_RESOURCES:
        downward_lf_reserves = model.Provide_LF_Downward_Reserve_From_Charge_MW[resource, timepoint]
    else:
        downward_lf_reserves = 0


    return downward_reg + downward_lf_reserves + model.Increase_Load_MW[resource, timepoint] \
           <= model.Max_Increase_Load_In_Timepoint_MW[resource, timepoint]


def intertimepoint_ramp_rule_lb(model, resource, timepoint, ramp_duration):
    """
    Inter-timepoint ramp constraints enforced for each resource in every timepoint.
    """

    if model.ramp_down_intertimepoint_limit_fraction[resource, ramp_duration]=="undefined":
        return pyo.Constraint.Skip
    else:

        if model.hour_of_day[timepoint] - ramp_duration < 0:
            timepoint_shift = model.timepoints_per_day
        else:
            timepoint_shift = 0

        net_power_provided = (model.Provide_Power_MW[resource, timepoint]
                                  - model.Provide_Power_MW[resource, timepoint - ramp_duration + timepoint_shift])

        return - (model.Operational_Capacity_MW[resource, model.period[timepoint]] *
                  model.ramp_down_intertimepoint_limit_fraction[resource, ramp_duration]) \
               <= net_power_provided


def intertimepoint_ramp_rule_ub(model, resource, timepoint, ramp_duration):
    """
    Inter-timepoint ramp constraints enforced for each resource in every timepoint.
    """

    if model.ramp_up_intertimepoint_limit_fraction[resource, ramp_duration] == "undefined":
        return pyo.Constraint.Skip
    else:

        if model.hour_of_day[timepoint] - ramp_duration < 0:
            timepoint_shift = model.timepoints_per_day
        else:
            timepoint_shift = 0

        net_power_provided = (model.Provide_Power_MW[resource, timepoint]
                                  - model.Provide_Power_MW[resource, timepoint - ramp_duration + timepoint_shift])

        return net_power_provided <= \
               (model.Operational_Capacity_MW[resource, model.period[timepoint]] *
                model.ramp_up_intertimepoint_limit_fraction[resource, ramp_duration])


def energy_shift_increase_load_adjacency_rule(model, resource, timepoint):
    """ resources that can shift load (storage, shift DR, etc.) need to make sure that
    Increase_Load_MW for a certain hour is no more than total Provide_Power_MW of adjacent hours as specified.
    
    Note this constraint is inclusive of the current timepoint.

    Args:
        model:
        resource that increase load & provide power:
        timepoint:

    Returns:
        Increase_load_mw <= sum(provide_power_mw) for adjacent hours
    """
    if model.adjacency_h[resource]==9999:
        return pyo.Constraint.Skip
    else:
        if resource in model.RESOURCES_WITH_ENERGY_IN_STORAGE:
            if (model.charging_efficiency[resource] != 1) or (model.discharging_efficiency[resource] != 1):
                raise ValueError("Charging and discharging efficiencies for resources "
                                 "with adjacency and energy in storage limits must equal 1.")
        rolling_sum_power = model.Provide_Power_MW[resource, timepoint]
        left_pointer = timepoint
        right_pointer = timepoint
        for step in range(0,model.adjacency_h[resource]):
            left_pointer = model.previous_timepoint[left_pointer]
            right_pointer = model.next_timepoint[right_pointer]
            rolling_sum_power += model.Provide_Power_MW[resource, left_pointer] \
                                 + model.Provide_Power_MW[resource, right_pointer]

        return model.Increase_Load_MW[resource, timepoint] <= rolling_sum_power


def energy_shift_provide_power_adjacency_rule(model, resource, timepoint):
    """ resources that can shift load (storage, shift DR, etc.) need to make sure that
        Provide_Power_MW for a certain hour is no more than total Increase_Load_MW of adjacent hours as specified.

        Args:
            model:
            resource that increase load & provide power:
            timepoint:

        Returns:
            Provide_Power_MW <= sum(Increase_Load_MW) for adjacent hours
        """
    if model.adjacency_h[resource] == 9999:
        return pyo.Constraint.Skip
    else:
        if resource in model.RESOURCES_WITH_ENERGY_IN_STORAGE:
            if (model.charging_efficiency[resource] != 1) or (model.discharging_efficiency[resource] != 1):
                raise ValueError("Charging and discharging efficiencies for resources "
                                 "with adjacency and energy in storage limits must equal 1.")
        rolling_sum_load_incr = model.Increase_Load_MW[resource, timepoint]
        left_pointer = timepoint
        right_pointer = timepoint
        for step in range(0, model.adjacency_h[resource]):
            left_pointer = model.previous_timepoint[left_pointer]
            right_pointer = model.next_timepoint[right_pointer]
            rolling_sum_load_incr += model.Increase_Load_MW[resource, left_pointer] \
                                 + model.Increase_Load_MW[resource, right_pointer]

        return model.Provide_Power_MW[resource, timepoint] <= rolling_sum_load_incr

# ### Storage ### #


def _enable_hourly_energy_tracking_validate(model, value, resource):
    """Hourly energy budget accounting is required for `STORAGE_RESOURCES` and optional for `RESOURCES_THAT_INCREASE_LOAD`. 
    
    Other resources must use the generic energy budget accounting constraint `Main_Energy_Budget_Constraint`
    """
    if model.storage[resource]:
        return value==True # must be true for storage resource
    elif not model.increases_load[resource]:
        return value==False # must be false if resource cannot increase load
    else:
        return True


def _load_only_validate(model, value, resource):
# storage, unit commitment, variable, thermal and generate_at_max resources cannot be load_only
    if (
            (model.storage[resource]) |
            (model.unit_commitment[resource]) |
            (model.variable[resource]) |
            (model.thermal[resource]) |
            (model.generate_at_max[resource])
    ):
        return value == False
    else:
        return True


def _increases_load_validate(model, value, resource):
# storage resources must be able to increase load
    if model.storage[resource]:
        return value == True  # must be True for storage resource
    else:
        return True

def min_energy_validation_rule(model, value, resource, timepoint):
    """
    Energy_In_Storage_MWh for storage resources should never be negative;
    Energy_In_Storage_MWh for load only resources should never be negative;
    Energy_In_Storage_MWh for flexible shift load resources can be negative
    """
    if resource in model.STORAGE_RESOURCES | model.LOAD_ONLY_RESOURCES:
        return value>=0.0
    else:
        return True


def available_energy_capacity_mwh_rule(model, resource, timepoint):
    if resource in model.STORAGE_RESOURCES:
        available_mwh = model.Total_Storage_Energy_Capacity_MWh[resource, model.period[timepoint]] \
           * model.available_capacity_fraction[resource, timepoint]
    else:
        if model.fixed_duration_h[resource] == 9999:
            available_mwh = float("inf")
        else:
            available_mwh = model.Available_Capacity_In_Timepoint_MW[resource, timepoint] \
                        * model.fixed_duration_h[resource]

    return available_mwh


def storage_energy_upper_rule(model, resource, timepoint):
    """No more total energy can be stored at any point that the total storage energy capacity.

    Args:
      model:
      timepoint:
      resource:

    Returns:

    """
    if resource not in model.STORAGE_RESOURCES and model.fixed_duration_h[resource] == 9999:
        return pyo.Constraint.Skip
    else:
        return model.Energy_in_Storage_MWh[resource, timepoint] \
            <= model.Available_Energy_Capacity_MWh[resource, timepoint]


def storage_energy_lower_rule(model, resource, timepoint):
    """
    Constrain minimum energy in storage in each timepoint (controlled by flexible param `min_energy_in_storage_in_timepoint_duration_h`).
    
    Note: maximum energy in storage is controlled via `Available_Energy_Capacity_MWh`, which is affected by either `available_capacity_fraction` for `STORAGE_RESOURCES` or `fixed_duration_h` all other `RESOURCES_WITH_ENERGY_IN_STORAGE`.
    """
    if model.min_energy_in_storage_in_timepoint_duration_h[resource, timepoint] == -9999:
        return pyo.Constraint.Skip
    else:
        return model.Energy_in_Storage_MWh[resource, timepoint] \
               >= model.Available_Capacity_In_Timepoint_MW[resource, timepoint] \
               * model.min_energy_in_storage_in_timepoint_duration_h[resource, timepoint]


def storage_energy_tracking_rule(model, storage_resource, timepoint):
    """The total energy in storage at the start of the next timepoint must equal
    the energy in storage at the start of the current timepoint
    plus charging that happened in the current timepoint, adjusted for the charging efficiency,
    minus discharging in the current timepoint, adjusted for the discharging efficiency.
    If providing reserve, assume a fraction of the time the reserves are dispatched, so add the respective additional
    energy stored from downward reserves (multiplied by the charging efficiency, as less is actually stored than
    needed from the grid) and additional energy discharged from upward reserves (divided by the discharging efficiency,
    as less is released to the grid than was stored).
    Assume providing frequency response or spinning reserve -- needed in rare contingency events
    -- will not affect the state of charge.
    Note that another constraint (storage_discharge_reserves_energy_rule) ensures that there is enough energy in the
    energy reservoir to dispatch spinning reserve for an entire timestep.
    energy_taken_offline_mwh & energy_taken_offline_mwh default to zero for most resources - these params are for
    resources that connect or disconnect from the grid or have end-use needs that require energy to be consumed or
    added to the grid at a certain time e.g. EVs connecting and disconnecting from grid

    Args:
      model:
      timepoint:
      storage_resource:

    Returns:

    """
    charging_mwh = 0.0
    discharging_mwh = 0.0
    energy_budget_current_tmp = 0.0
    daily_energy_budget_change_mwh = 0.0

    # Charging / Increase Load
    charging_mwh += model.Increase_Load_MW[storage_resource, timepoint]
    # Reserves from charging
    if storage_resource in model.REGULATION_RESERVE_RESOURCES:
        charging_mwh += (model.Provide_Downward_Reg_From_Charge_MW[storage_resource, timepoint]
                        - model.Provide_Upward_Reg_From_Charge_MW[storage_resource, timepoint]) \
            * model.reg_dispatch_fraction        
    if storage_resource in model.LOAD_FOLLOWING_RESERVE_RESOURCES:
        charging_mwh += (model.Provide_LF_Downward_Reserve_From_Charge_MW[storage_resource, timepoint]
                        - model.Provide_LF_Upward_Reserve_From_Charge_MW[storage_resource, timepoint]) \
            * model.lf_reserve_dispatch_fraction
        
    # Discharging / Provide Power
    if storage_resource in model.RESOURCES_WITH_PROVIDE_POWER:
        discharging_mwh += model.Provide_Power_MW[storage_resource, timepoint]
        # Reserves from discharging
        if storage_resource in model.REGULATION_RESERVE_RESOURCES:
            discharging_mwh += (model.Provide_Upward_Reg_MW[storage_resource, timepoint]
                                - model.Provide_Downward_Reg_MW[storage_resource, timepoint]) \
                * model.reg_dispatch_fraction      
        if storage_resource in model.LOAD_FOLLOWING_RESERVE_RESOURCES:
            discharging_mwh += (model.Provide_LF_Upward_Reserve_MW[storage_resource, timepoint]
                                - model.Provide_LF_Downward_Reserve_MW[storage_resource, timepoint]) \
                * model.lf_reserve_dispatch_fraction
            
        if model.energy_budget_cap_factor[storage_resource, model.day[timepoint]] != 'undefined':
            energy_budget_current_tmp = model.Operational_Capacity_MW[storage_resource, model.period[timepoint]] \
                                        * model.energy_budget_cap_factor[storage_resource, model.day[timepoint]]
        if model.interday_energy_sharing:
            if storage_resource in model.RESOURCES_WITH_MULTIDAY_SHARING:
                daily_energy_budget_change_mwh = \
                    model.Daily_Energy_Budget_Change_MWh[storage_resource, model.period[timepoint], model.day[timepoint]]

    return model.Energy_in_Storage_MWh[storage_resource, model.next_timepoint[timepoint]] \
        == model.Energy_in_Storage_MWh[storage_resource, timepoint] \
        + charging_mwh * model.charging_efficiency[storage_resource] \
        + daily_energy_budget_change_mwh \
            / model.timepoints_per_day \
        + energy_budget_current_tmp \
        - discharging_mwh / model.discharging_efficiency[storage_resource] \
        - model.energy_taken_offline_mwh[storage_resource, timepoint] \
        + model.energy_brought_online_mwh[storage_resource, timepoint]


def storage_discharge_reserves_energy_rule(model, storage_resource, timepoint):
    """If called upon to provide upward reserve, storage must be have enough energy available to sustain output for
    the duration of the timepoint. Because of the discharge (in)efficiency, it will actually release less energy to the
    grid than is actually available in storage, so multiply by the discharge efficiency.
    Load-following, regulation, and spin are assumed to need to be available for the full duration of the timepoint
    if called upon.
    Ignore frequency response provision on the assumption that it does not significantly affect state of charge.

    Args:
      model:
      timepoint:
      storage_resource:

    Returns:

    This rule covers only the reserves provided by discharge.

    """
    if (
        storage_resource in model.LOAD_FOLLOWING_RESERVE_RESOURCES or
        storage_resource in model.REGULATION_RESERVE_RESOURCES or storage_resource in model.SPINNING_RESERVE_RESOURCES
    ):

        if storage_resource in model.REGULATION_RESERVE_RESOURCES:
            upward_reg = model.Provide_Upward_Reg_MW[storage_resource, timepoint]
        else:
            upward_reg = 0

        if storage_resource in model.LOAD_FOLLOWING_RESERVE_RESOURCES:
            upward_lf_reserves = model.Provide_LF_Upward_Reserve_MW[storage_resource, timepoint]
        else:
            upward_lf_reserves = 0

        if storage_resource in model.SPINNING_RESERVE_RESOURCES:
            spin = model.Provide_Spin_MW[storage_resource, timepoint]
        else:
            spin = 0

        return model.Provide_Power_MW[storage_resource, timepoint] \
            + upward_reg + upward_lf_reserves + spin \
            <= model.Energy_in_Storage_MWh[storage_resource, timepoint] \
            * model.discharging_efficiency[storage_resource]
    else:
        return pyo.Constraint.Skip


def storage_charge_reserves_energy_rule(model, storage_resource, timepoint):
    """If called upon to provide downward reserve, storage must be have enough energy storage capacity to store the
    additional energy from the new set point. Because of the charging (in)efficiency, it will have to get more energy
    from the grid than it will be able to put in storage, so divide by the charging efficiency.
    Load-following and regulation are assumed to have to be available for the full duration of the timepoint if called
    upon.

    Args:
      model:
      timepoint:
      storage_resource:

    Returns:

    This rule covers only the reserves provided by charge.

    """
    if (
        storage_resource in model.LOAD_FOLLOWING_RESERVE_RESOURCES or
        storage_resource in model.REGULATION_RESERVE_RESOURCES
    ):

        if storage_resource in model.REGULATION_RESERVE_RESOURCES:
            downward_reg = model.Provide_Downward_Reg_From_Charge_MW[storage_resource, timepoint]
        else:
            downward_reg = 0

        if storage_resource in model.LOAD_FOLLOWING_RESERVE_RESOURCES:
            downward_lf_reserves = model.Provide_LF_Downward_Reserve_From_Charge_MW[storage_resource, timepoint]
        else:
            downward_lf_reserves = 0

        return model.Increase_Load_MW[storage_resource, timepoint] \
            + downward_reg + downward_lf_reserves \
            <= (model.Available_Energy_Capacity_MWh[storage_resource, timepoint] -
                model.Energy_in_Storage_MWh[storage_resource, timepoint]) \
            / model.charging_efficiency[storage_resource]
    else:
        return pyo.Constraint.Skip






def paired_battery_charging(model, r, t):
    """Constrain a paired storage resource by the energy production of its paired generator.

    Args:
      model:
      r:
      t:

    Returns:

    """
    storage_charging = model.Increase_Load_MW[r, t]
    if r in model.REGULATION_RESERVE_RESOURCES:
        storage_charging += model.Provide_Downward_Reg_From_Charge_MW[r, t] * model.reg_dispatch_fraction

    if r in model.LOAD_FOLLOWING_RESERVE_RESOURCES:
        storage_charging += model.Provide_LF_Downward_Reserve_From_Charge_MW[r, t] * model.lf_reserve_dispatch_fraction

    return storage_charging <= model.Provide_Power_MW[model.paired_supply_resource[r], t]


def discharge_efficiency_validate(model,discharge_efficiency,resource):
    """
    Discharge efficiency is applied when a resource provides power.
    Load only resources never provide power, so force efficieny to be 100% - a stand-in for undefined.
    Resources with energy budgets must have discharge efficiency == 100% for RPS accounting.
    """
    
    if resource in model.LOAD_ONLY_RESOURCES:
        return discharge_efficiency == 1.0
    elif all(
        model.energy_budget_cap_factor[resource,day] == 'undefined'
        for day in model.DAYS
    ):
        return 0.0 <= discharge_efficiency <= 1.0
    else:
        return discharge_efficiency == 1.0


def _can_produce_less_than_energy_budget_validate(model,value,resource):
    """
    Storage resources should not be allowed to produce less than energy budget
    """
    if resource in model.STORAGE_RESOURCES:
        return value == False #must be false if storage
    else:
        return True


def _energy_storage_charge_capacity_cost_validate(model,value,resource,index):
    """
    Resources with a nonzero charge capacity must have an associated cost (to prevent forgetting to enter this)

    More technical definition: if the cost is zero (the default value), then the planned storage charge capacity must be
    inf (the default value)
    """
    if value == 0:
        return all(model.planned_storage_charge_capacity_mw[resource,period]==float("inf") for period in model.PERIODS)
    else:
        return True

# ##### Power Balance and System Operational Constraints ##### #


def zonal_power_balance_rule(model, zone, timepoint):
    """The sum of all in-zone power production and net transmission flow equals the zone's load in each timepoint.
    Resources may increase load; the demand from these resources is added to the zone's load in each timepoint.
    Storage resources both increase load and provide power, and are therefore included on both sides of the constraint.
    Two slack variables for unserved energy and overgeneration are included.
    Scheduled renewable curtailment is not explicitly included here because Provide_Power_MW
    will be less than the resource's potential power production during times of scheduled curtailment

    Args:
      model:
      timepoint:
      zone:

    Returns:

    """

    # energy efficiency
    ee_load_reduction = float()
    for resource in model.RESOURCES:
        if model.zone[resource] == zone:
            if resource in model.EE_PROGRAMS:
                ee_load_reduction += model.EE_Reduced_Load_FTM_MW[resource, timepoint]

    # Imports/exports
    # Transmit_Power_Unspecified is used here because dedicated imports
    # are modeled as supplying power to the "to" zone directly,
    # so power from dedicated import resources will be included with resources above
    imports_exports = 0
    for line in model.TRANSMISSION_LINES:
        if model.transmission_to[line] == zone or model.transmission_from[line] == zone:
            if model.transmission_to[line] == zone:
                imports_exports += model.Transmit_Power_Unspecified_MW[line, timepoint]
            elif model.transmission_from[line] == zone:
                imports_exports -= model.Transmit_Power_Unspecified_MW[line, timepoint]

    # Imports and exports through semi storage zones
    # The positive direction for SSZ_Transmit_Power_MW is from the zone to storage semi zones,
    # so the direction is negative
    semi_storage_zones_imports_exports = float()
    if model.allow_semi_storage_zones:
        for ssz in model.SEMI_STORAGE_ZONES:
            if model.ssz_from_zone[ssz] == zone:
                semi_storage_zones_imports_exports -= model.SSZ_Transmit_Power_MW[ssz, timepoint]

    if model.allow_unserved_energy:
        unserved_energy = model.Unserved_Energy_MW[zone, timepoint]
    else:
        unserved_energy = 0.0

    # provide power and increase load are already summed for each zone, so use the expressions here.
    # Unserved energy & overgeneration for infeasibility diagnosis
    return model.Zonal_Provide_Power_MW[zone, timepoint] \
           + imports_exports \
           + semi_storage_zones_imports_exports \
           - model.Overgeneration_MW[zone, timepoint] \
           == model.input_load_mw[zone, timepoint] \
              + model.Zonal_Increase_Load_MW[zone, timepoint] \
              - ee_load_reduction \
              - unserved_energy


def transmission_min_flow_rule(model, line, timepoint):
    """Transmission flows must obey flow limits on each line.

    Args:
      model:
      timepoint:
      line:

    Returns:

    """
    # Default tx bound
    min_flow_bound = model.min_flow_planned_mw[line]
    # If allow_tx_build, incremental capacity is added onto existing
    if model.allow_tx_build and line in model.TRANSMISSION_LINES_NEW:
        min_flow_bound -= model.New_Tx_Total_Installed_Capacity_MW[line, model.period[timepoint]]

    return min_flow_bound <= model.Transmit_Power_MW[line, timepoint]


def transmission_max_flow_rule(model, line, timepoint):
    """Transmission flows must obey flow limits on each line.

    Args:
      model:
      timepoint:
      line:

    Returns:

    """
    # Default tx bound
    max_flow_bound = model.max_flow_planned_mw[line]
    # If allow_tx_build, incremental capacity is added onto existing
    if model.allow_tx_build and line in model.TRANSMISSION_LINES_NEW:
        max_flow_bound += model.New_Tx_Total_Installed_Capacity_MW[line, model.period[timepoint]]

    return model.Transmit_Power_MW[line, timepoint] <= max_flow_bound


def transmission_ramp_down_rule(model, line, tmp, ramp_duration):
    """Multi-hour intertie ramp down constraints enforced in every timepoint.

    Args:
      model:
      tmp:
      line:
      ramp_duration:

    Returns:

    """
    if model.hour_of_day[tmp] - ramp_duration < 0:
        timepoint_shift = model.timepoints_per_day
    else:
        timepoint_shift = 0

    # Default tx bound
    max_flow_bound = (model.max_flow_planned_mw[line] - model.min_flow_planned_mw[line])
    # If allow_tx_build, incremental capacity is added onto existing
    if model.allow_tx_build and line in model.TRANSMISSION_LINES_NEW:
        # add the incremental capacity in both the positive and negative directions
        max_flow_bound += 2 * model.New_Tx_Total_Installed_Capacity_MW[line, model.period[tmp]]

    return (
        -model.flow_ramp_down_limit_fraction[line,
                                             ramp_duration] * max_flow_bound <=
        model.Transmit_Power_MW[line,
                                tmp] - model.Transmit_Power_MW[line,
                                                               tmp - ramp_duration + timepoint_shift]
    )


def transmission_ramp_up_rule(model, line, tmp, ramp_duration):
    """Multi-hour intertie ramp up constraints enforced in every timepoint.

    Args:
      model:
      tmp:
      line:
      ramp_duration:

    Returns:

    """
    if model.hour_of_day[tmp] - ramp_duration < 0:
        timepoint_shift = model.timepoints_per_day
    else:
        timepoint_shift = 0

    # Default tx bound
    max_flow_bound = (model.max_flow_planned_mw[line] - model.min_flow_planned_mw[line])
    # If allow_tx_build, incremental capacity is added onto existing
    if model.allow_tx_build and line in model.TRANSMISSION_LINES_NEW:
        max_flow_bound += 2 * model.New_Tx_Total_Installed_Capacity_MW[line, model.period[tmp]]

    return (
        model.Transmit_Power_MW[line,
                                tmp] - model.Transmit_Power_MW[line,
                                                               tmp - ramp_duration + timepoint_shift] <=
        model.flow_ramp_up_limit_fraction[line,
                                          ramp_duration] * max_flow_bound
    )


def simultaneous_transmission_flows_rule(model, group, timepoint):
    """The sum of flows on a group of lines (in a certain direction) cannot exceed a pre-specified limit.

    Args:
      model:
      timepoint:
      group:

    Returns:

    """
    sim_flow = float()
    for (g, line) in model.SIMULTANEOUS_FLOW_GROUP_LINES:
        if g == group:
            sim_flow += model.Transmit_Power_MW[line, timepoint] * model.direction[g, line]
    return sim_flow <= model.simultaneous_flow_limit_mw[group, model.period[timepoint]]


# ### Operational Reserves ### #
def spinning_reserve_req(model, timepoint):
    """

    Args:
      model:
      timepoint:

    Returns:

    """
    spinning_reserve_requirement = float()
    for zone in model.ZONES:
        # how much does energy efficiency reduce demand?
        # Will be used to adjust the spinning reserve requirement
        ee_load_reduction = 0
        for resource in model.EE_PROGRAMS:
            if model.zone[resource] == zone:
                ee_load_reduction += model.EE_Reduced_Load_FTM_MW[resource, timepoint]

        spinning_reserve_requirement += model.spin_reserve_fraction_of_load[zone] * \
            (model.input_load_mw[zone, timepoint] - ee_load_reduction)

    # add spin reserved specified by MW (will only be nonzero if spin_reserve_fraction_of_load is zero, by validation)
    spinning_reserve_requirement += model.spin_reserve_req_mw[timepoint]

    return spinning_reserve_requirement


def spinning_reserve_rule(model, timepoint):
    """Meet the spinning reserve requirement in each timepoint.

    Args:
      model:
      timepoint:

    Returns:

    """
    spinning_reserve_provision = float()
    for resource in model.SPINNING_RESERVE_RESOURCES:
        spinning_reserve_provision += model.Provide_Spin_Net_MW[resource, timepoint]

    return (spinning_reserve_provision + model.Spin_Violation_MW[timepoint] == model.Spinning_Reserve_Req_MW[timepoint])


def upward_regulation_rule(model, timepoint):
    """Meet the regulation requirement in each timepoint.

    Args:
      model:
      timepoint:

    Returns:
      rule ensuring upward regulation requirement is met

    """
    upward_regulation_provision = float()
    for resource in model.REGULATION_RESERVE_RESOURCES:
        upward_regulation_provision += model.Provide_Upward_Reg_Net_MW[resource, timepoint]

    return upward_regulation_provision + model.Upward_Reg_Violation_MW[timepoint] \
        == model.upward_reg_req[timepoint]


def downward_regulation_rule(model, timepoint):
    """Meet the downward regulation requirement in each timepoint.
    Currently variable renewables are assumed to not provide regulation.

    Args:
      model:
      timepoint:

    Returns:

    """
    downward_regulation_provision = float()
    for resource in model.REGULATION_RESERVE_RESOURCES:
        downward_regulation_provision += model.Provide_Downward_Reg_Net_MW[resource, timepoint]

    return downward_regulation_provision + model.Downward_Reg_Violation_MW[timepoint] \
        == model.downward_reg_req[timepoint]


def upward_combined_lf_reserve_req(model, timepoint):
    """Add endogenous LF requirement to baseline input LF requirement if applicable.

    The expression assumes all variable resources within the set of load following zones contribute to the endogenous LF requirement
    and scales the incremental system load following need linearly with Operational_Capacity_MW

    Args:
      model:
      timepoint:

    Returns:

    """
    req = model.upward_lf_reserve_req[timepoint]

    for resource in model.VARIABLE_RESOURCES:
        if model.zone[resource] in model.LOAD_FOLLOWING_ZONES:
            req += (
                model.Operational_Capacity_MW[resource,
                                              model.period[timepoint]] *
                model.resource_upward_lf_req[resource,
                                             model.day[timepoint],
                                             model.hour_of_day[timepoint]]
            )
    return req


def downward_combined_lf_reserve_req(model, timepoint):
    """Add endogenous LF requirement to baseline input LF requirement if applicable.

    The expression assumes all variable resources within the set of load following zones contribute to the endogenous LF requirement
    and scales the incremental system load following need linearly with Operational_Capacity_MW

    Args:
      model:
      timepoint:

    Returns:

    """
    req = model.downward_lf_reserve_req[timepoint]

    for resource in model.VARIABLE_RESOURCES:
        if model.zone[resource] in model.LOAD_FOLLOWING_ZONES:
            req += (
                model.Operational_Capacity_MW[resource,
                                              model.period[timepoint]] *
                model.resource_downward_lf_req[resource,
                                               model.day[timepoint],
                                               model.hour_of_day[timepoint]]
            )
    return req


def downward_load_following_reserve_rule(model, timepoint):
    """Meet downward load following reserve requirement by incurring subhourly curtailment.

    Downward LF not provided by firm resources or taken as a violation must be provided by curtailable
    variable resources.

    Args:
      model:
      timepoint:

    Returns:

    """
    firm_downward_lf_provision_mw = sum(
        model.Provide_LF_Downward_Reserve_Net_MW[r,
                                             timepoint] for r in model.LOAD_FOLLOWING_RESERVE_RESOURCES
    )
    # make a separate Var for Sub_Curtailment / lf_dispatch
    return (
        model.Variable_Resource_Provide_Downward_LF_MW[timepoint] + model.Downward_LF_Reserve_Violation_MW[timepoint] +
        firm_downward_lf_provision_mw == model.Downward_Load_Following_Reserve_Req[timepoint]
    )


def meet_upward_load_following_reserve_rule(model, timepoint):
    """Ensure that upward load following reserves are met.

    Args:
      model:
      timepoint:

    Returns:

    """
    firm_upward_lf_provision_mw = sum(
        model.Provide_LF_Upward_Reserve_Net_MW[r,
                                           timepoint] for r in model.LOAD_FOLLOWING_RESERVE_RESOURCES
    )

    return (
        firm_upward_lf_provision_mw + model.Variable_Resource_Provide_Upward_LF_MW[timepoint] +
        model.Upward_LF_Reserve_Violation_MW[timepoint] == model.Upward_Load_Following_Reserve_Req[timepoint]
    )


def upward_load_following_reserve_rule(model, timepoint):
    """Limit the available upward LF from variable resources as a fraction of total scheduled curtailment.

    Variable resources can provide as much upward LF as they've scheduled for curtailment,
    even though only a fraction of the bid is dispatched (e.g., 20%)

    Args:
      model:
      timepoint:

    Returns:

    """
    if model.variable_resources_upward_lf:
        curtailment_mw = sum(
            model.Scheduled_Curtailment_MW[r,
                                           timepoint] for r in model.LOAD_FOLLOWING_ZONE_CURTAILABLE_RESOURCES
        )
    else:
        curtailment_mw = 0

    return (
        model.Variable_Resource_Provide_Upward_LF_MW[timepoint] <=
        model.var_rnw_available_for_lf_reserves * curtailment_mw
    )


def variable_rnw_down_reserve_availability_rule(model, timepoint):
    """Limit the available downward LF from variable resources as a fraction of total variable resource power production.

    Args:
      model:
      timepoint:

    Returns:

    """
    total_variable_renewables = float()
    for resource in model.LOAD_FOLLOWING_ZONE_CURTAILABLE_RESOURCES:
        total_variable_renewables += model.Provide_Power_MW[resource, timepoint]

    return (
        model.Variable_Resource_Provide_Downward_LF_MW[timepoint] <=
        model.var_rnw_available_for_lf_reserves * total_variable_renewables
    )


def max_upward_lf_from_variable_rule(model, timepoint):
    """Limit the fraction of upward LF reserves that can be provided by variable resources.

    Args:
      model:
      timepoint:

    Returns:

    """
    return (
        model.Variable_Resource_Provide_Upward_LF_MW[timepoint] <=
        model.max_var_rnw_lf_reserves * model.Upward_Load_Following_Reserve_Req[timepoint]
    )


def variable_rnw_downward_lf_reserve_limit_rule(model, timepoint):
    """Limit the fraction of the reserve requirement can be met with variable generation.

    Args:
      model:
      timepoint:

    Returns:

    """
    return (
        model.Variable_Resource_Provide_Downward_LF_MW[timepoint] <=
        model.max_var_rnw_lf_reserves * model.Downward_Load_Following_Reserve_Req[timepoint]
    )


def total_frequency_response_rule(model, timepoint):
    """Available headroom from resources contributing to the total frequency response requirement
    must be greater than the requirement.
    In other operational constraints, frequency response is modeled as additive to upward reserves,
    which means that provision of the two services is not allowed to be overlapping.
    Skip constraint if requirement is 0.

    Args:
      model:
      timepoint:

    Returns:

    """
    if model.freq_resp_total_req_mw[timepoint] == 0:
        return pyo.Constraint.Skip
    else:
        frequency_response_provision = float()
        for resource in model.TOTAL_FREQ_RESP_RESOURCES:
            frequency_response_provision += model.Provide_Frequency_Response_Net_MW[resource, timepoint]

        return frequency_response_provision >= model.freq_resp_total_req_mw[timepoint]


def partial_frequency_response_rule(model, timepoint):
    """Available headroom from resources contributing to the partial frequency response requirement
    must be greater than the requirement.
    Partial frequency response headroom can be shared with total frequency response headroom,
    but not with other upward reserves.

    Args:
      model:
      timepoint:

    Returns:

    """

    if model.freq_resp_partial_req_mw[timepoint] == 0:
        return pyo.Constraint.Skip
    else:
        frequency_response_provision = float()
        for resource in model.PARTIAL_FREQ_RESP_RESOURCES:
            frequency_response_provision += model.Provide_Frequency_Response_Net_MW[resource, timepoint]

        return frequency_response_provision >= model.freq_resp_partial_req_mw[timepoint]


def minimum_local_committed_generation_rule(model, timepoint):
    """A certain amount of eligible generation must be committed in all hours
    from resources in the set of MINIMUM_GENERATION_RESOURCES
    for any timepoint that has min_gen_committed_mw > 0.

    Args:
      model:
      timepoint:

    Returns:

    """
    if model.min_gen_committed_mw[timepoint] == 0:
        return pyo.Constraint.Skip
    else:
        local_generation = float()
        for resource in model.MINIMUM_GENERATION_RESOURCES:
            local_generation += model.Commit_Capacity_MW[resource, timepoint]

        return local_generation >= model.min_gen_committed_mw[timepoint]


# ### Deliverability ### #
def deliverability_def_rule(model, resource, period):
    """Split newly installed resource capacity between delivered and undelivered capacity
    Deliverability functionality currently includes only resources that
    A) Do not have planned capacity (planned_installed_capacity_mw = 0 for all periods)
    B) Cannot retire (not in CAN_RETIRE_RESOURCES)
    C) Can be newly built (in NEW_BUILD_RESOURCES)
    These restrictions are enforced via the validation rule "deliverability_resources_validation_rule"
    Cumulative_New_Installed_Capacity_MW is used in this constraint because
    it is consistent with the three conditions above.
    Note that as currently defined, Operational_New_Capacity_MW = Cumulative_New_Installed_Capacity_MW
    for any deliverability resource

    Args:
      model:
      period:
      resource:

    Returns:

    """
    return model.Delivered_Capacity_MW[resource, period] \
        + model.Undelivered_Capacity_MW[resource, period] \
        == model.Cumulative_New_Installed_Capacity_MW[resource, period]


def delivered_capacity_period_build_rule(model, resource, period):
    """Delivered_Capacity_MW cannot decrease over time.
    An identical constraint exists for Undelivered_Capacity_MW

    Args:
      model:
      period:
      resource:

    Returns:

    """
    if period == model.first_period:
        return pyo.Constraint.Skip
    else:
        return model.Delivered_Capacity_MW[resource, period] \
            >= model.Delivered_Capacity_MW[resource, find_prev_period(model, period)]


def undelivered_capacity_period_build_rule(model, resource, period):
    """Undelivered_Capacity_MW cannot decrease over time.
    An identical constraint exists for Delivered_Capacity_MW

    Args:
      model:
      period:
      resource:

    Returns:

    """
    if period == model.first_period:
        return pyo.Constraint.Skip
    else:
        return model.Undelivered_Capacity_MW[resource, period] \
            >= model.Undelivered_Capacity_MW[resource, find_prev_period(model, period)]

def new_deliverability_capacity_tracking_rule(model, deliverability_zone, period):
    """Track how much new deliverability capacity we have in each deliverability zone in each period.

    Args:
      model:
      deliverability_zone:
      period:

    Returns:

    """
    new_deliverability_capacity = float()
    for vintage in model.VINTAGES:
        if vintage <= period:
            new_deliverability_capacity += model.New_Deliverability_Incremental_MW[deliverability_zone, vintage]

    return new_deliverability_capacity

def deliverability_cost_annual_in_period_dollars_rule(model, deliverability_zone, period):
    """Add levelized capital costs for all deliverability_zone additions from previous periods and the current period

    Args:
      model:
      deliverability_zone:
      period:

    Returns:

    """
    deliverability_cost_total = 0.0

    for vintage in model.VINTAGES:
        if period >= vintage:
            deliverability_cost_total += model.New_Deliverability_Incremental_MW[deliverability_zone, vintage] \
                                  * model.new_deliverability_cost_per_mw_yr[deliverability_zone, vintage]

    return deliverability_cost_total



# ############ Transmission Direction #############
# The following two constraints define the positive and negative direction versions of Transmit_Power_MW,
# which can then be multiplied by directional hurdle, emissions import rate, or other value.
# Important note: any value multiplied by these variables cannot provide an incentive
# for the directional variable to exceed the absolute value of the non-directional variable in the same direction.
# This means that, for example, a negative cost hurdle rate on the positive direction would give an incentive for
# Transmit_Power_Positive_Direction_MW to be much larger than Transmit_Power_MW
# which is not acceptable because the variables are supposed to be equal to each other when in the same direction.
def positive_direction_transmit_power_rule(model, line, timepoint):
    """

    Args:
      model:
      timepoint:
      line:

    Returns:

    """
    return model.Transmit_Power_Unspecified_Positive_Direction_MW[line, timepoint] \
        >= model.Transmit_Power_Unspecified_MW[line, timepoint]


def negative_direction_transmit_power_rule(model, line, timepoint):
    """

    Args:
      model:
      timepoint:
      line:

    Returns:

    """
    return model.Transmit_Power_Unspecified_Negative_Direction_MW[line, timepoint] \
        >= -model.Transmit_Power_Unspecified_MW[line, timepoint]


# Conventional (shed) demand response
def conventional_dr_max_annual_availability_rule(model, dr_resource, period):
    """Total quantity of energy that can be interrupted (permanently shed) in each period.

    Args:
      model:
      period:
      dr_resource:

    Returns:

    """
    conventional_dr_dispatch = 0.0

    # Sum the MWh of conventional DR dispatch in the period
    for timepoint in model.TIMEPOINTS:
        if model.period[timepoint] == period:
            conventional_dr_dispatch += \
                model.Provide_Power_MW[dr_resource, timepoint] * model.day_weight[model.day[timepoint]]

    # Limit the MWh of conventional DR dispatch to the MW capacity
    # multiplied by the hours per year that the DR resource can be called (the RHS is MWh of availability)
    return conventional_dr_dispatch \
           <= model.Operational_Capacity_MW[dr_resource, period] \
           * model.conventional_dr_availability_hours_per_year[dr_resource, period]


def max_daily_provide_power_rule(model, resource, period, day):
    """Constrain total provide power within a day

    Args:
      model:
      dr_resource:
      day:
      period:

    Returns:

    """
    if model.max_daily_capacity_factor[resource, period, day] == 1.0:
        return pyo.Constraint.Skip
    else:
        return (
            sum(
                model.Provide_Power_MW[resource,timepoint]
                for timepoint in model.TIMEPOINTS
                if model.day[timepoint] == day and model.period[timepoint] == period
            ) <= model.Operational_Capacity_MW[resource,
                                               period] * model.max_daily_capacity_factor[resource, period, day] *
            model.timepoints_per_day
        )


def min_daily_increase_load_rule(model, resource, period, day):
    """Set a minimum level of increase load that must be met within each day.
    Enforced on the variable Increase_Load_MW, which is the gross amount of load increased.
    If a resource can also provide power, this constraint doesn't net off any power provided
    from the daily increase load minimum.

    Args:
      model:
      dr_resource:
      day:
      period:

    Returns:

    """
    if model.min_daily_increase_load_cap_factor[resource, period, day] == 0.0:
        return pyo.Constraint.Skip
    else:
        return (
            sum(
                model.Increase_Load_MW[resource,timepoint]
                for timepoint in model.TIMEPOINTS
                if model.day[timepoint] == day and model.period[timepoint] == period
            ) >= model.Operational_Capacity_MW[resource,
                                               period] * model.min_daily_increase_load_cap_factor[resource, period, day] *
            model.timepoints_per_day
        )


# ##### Constraints for multi-day energy sharing ##### #
def net_zero_energy_sharing_rule(model, resource, period, energy_sharing_interval):
    """
    Make sure the total energy within the sharing group stays the same before and after sharing
    """
    net_energy_moved = 0.0
    for day in model.DAYS:
        if model.interday_energy_sharing_interval_id[day] == energy_sharing_interval:
            net_energy_moved += (
                model.Daily_Energy_Budget_Change_MWh[resource, period, day] *
                model.day_weight[day]
            )

    return net_energy_moved == 0.0


def min_daily_changes_rule(model, resource, period, day):
    """
    Constrain daily energy budget changes to be within a user defined range
    """

    return (
        - (model.Operational_Capacity_MW[resource, period] *
           model.daily_max_energy_budget_decrease_hours[resource, day])
        <=
        model.Daily_Energy_Budget_Change_MWh[resource, period, day]
    )


def max_daily_changes_rule(model, resource, period, day):
    """
    Constrain daily energy budget changes to be within a user defined range
    """
    return (
        model.Daily_Energy_Budget_Change_MWh[resource, period, day]
        <=
        model.Operational_Capacity_MW[resource, period] *
        model.daily_max_energy_budget_increase_hours[resource, day]
    )


def define_absolute_energy_moved_rule(model, resource, period, day):
    """Get the only positive values of energy increase.

    Given the net zero constraint, net_zero_energy_sharing_rule, we know that exactly
    half of the energy will be increases and half will be decreases. In plain English,
    this constraint now means that you have x GWh of energy to shift to other parts
    of the year.
    """
    return (
        model.Positive_Energy_Budget_Moved_MWh[resource, period, day]
        >=
        model.Daily_Energy_Budget_Change_MWh[resource, period, day]
    )


def max_absolute_energy_moved_rule(model, resource, period, energy_sharing_interval):
    """Constraint the total amount of energy being moved
    """
    absolute_energy_moved = 0.0
    for day in model.DAYS:
        if model.interday_energy_sharing_interval_id[day] == energy_sharing_interval:
            absolute_energy_moved += (
                model.Positive_Energy_Budget_Moved_MWh[resource, period, day] *
                model.day_weight[day]
            )

    return absolute_energy_moved <= (
        model.Operational_Capacity_MW[resource, period] *
        model.max_energy_to_share_hours[resource, energy_sharing_interval]
    )


# initiate constraints


def maximum_ee_investment_in_period_rule(model, resource, period):
    """
    Limit EE investment in each period to max_investment_in_period_aMW
    """

    return model.Build_Capacity_MW[resource, period] <= model.max_investment_in_period_aMW[resource, period]


# ------------------------------------ Semi storage zones related constraints ----------------------------------------
def ssz_energy_net_0_rule(model, ssz_zone, period):
    """
    Make sure the energy transmit to and from ssz_zone within a period is net 0
    day weights are applied because the energy accounting.
    Emissions from semi storage zones are assumed to be zero,
    so semi storage zones don't appear anywhere in emissions constraints.
    """

    total_transmit_power_mw = 0.0
    for timepoint in model.TIMEPOINTS:
        if model.period[timepoint] == period:
            total_transmit_power_mw += model.SSZ_Transmit_Power_MW[ssz_zone, timepoint] \
                                       * model.day_weight[model.day[timepoint]]

    return total_transmit_power_mw == 0.0


# The following two constraints define the positive and negative direction versions of SSZ_Transmit_Power_MW,
# which can then be multiplied by directional hurdle, emissions import rate, hurdle rates, or other value.
# Important note: any value multiplied by these variables cannot provide an incentive
# for the directional variable to exceed the absolute value of the non-directional variable in the same direction.
# This means that, for example, a negative cost hurdle rate on the positive direction would give an incentive for
# Transmit_Power_Positive_Direction_MW to be much larger than Transmit_Power_MW
# which is not acceptable because the variables are supposed to be equal to each other when in the same direction.


def ssz_transmit_power_positive_definition_rule(model, ssz_zone, timepoint):
    return model.SSZ_Positive_Transmit_Power_MW[ssz_zone, timepoint] \
        >= model.SSZ_Transmit_Power_MW[ssz_zone, timepoint]


def ssz_transmit_power_negative_definition_rule(model, ssz_zone, timepoint):
    return model.SSZ_Negative_Transmit_Power_MW[ssz_zone, timepoint] \
        >= - model.SSZ_Transmit_Power_MW[ssz_zone, timepoint]


def ssz_tx_limit_max_rule(model, ssz_zone, timepoint):
    """
    Constraint the transmitting power is below the transmission limits
    """
    return model.SSZ_Transmit_Power_MW[ssz_zone, timepoint] <= model.ssz_max_flow_mw[ssz_zone, model.period[timepoint]]


def ssz_tx_limit_min_rule(model, ssz_zone, timepoint):
    """
    Constraint the transmitting power is below the transmission limits
    """
    return model.ssz_min_flow_mw[ssz_zone, model.period[timepoint]] \
        <= model.SSZ_Transmit_Power_MW[ssz_zone, timepoint]
