diff --git a/gridpath/system/water/water_flows.py b/gridpath/system/water/water_flows.py index 44c3ee2c1..a91dbb2e3 100644 --- a/gridpath/system/water/water_flows.py +++ b/gridpath/system/water/water_flows.py @@ -77,9 +77,10 @@ def water_link_departure_arrival_tmp_init(mod): return wl_dep_arr_tmp + m.TMPS_AND_OUTSIDE_HORIZON = Set(initialize=m.TMPS | {"tmp_outside_horizon"}) m.WATER_LINK_DEPARTURE_ARRIVAL_TMPS = Set( dimen=3, - within=m.WATER_LINKS * m.TMPS * m.TMPS, + within=m.WATER_LINKS * m.TMPS * m.TMPS_AND_OUTSIDE_HORIZON, initialize=water_link_departure_arrival_tmp_init, ) @@ -113,13 +114,14 @@ def max_flow_rule(mod, wl, dep_tmp, arr_tmp): def determine_arrival_timepoint(mod, tmp, travel_time_hours): # If this is the last timepoint of a linear horizon, there are no # timepoints to check + # TODO: check if this makes sense if check_if_boundary_type_and_last_timepoint( mod=mod, tmp=tmp, balancing_type=mod.water_system_balancing_type, boundary_type="linear", ): - tmp_to_check = None + tmp_to_check = "tmp_outside_horizon" elif check_if_boundary_type_and_last_timepoint( mod=mod, tmp=tmp, @@ -128,8 +130,12 @@ def determine_arrival_timepoint(mod, tmp, travel_time_hours): ): # TODO: add linked tmp_to_check = None - # Otherwise, we can start searching + # If travel time is less than the hours in the starting timepoint, we stay + # in the same timepoint + elif travel_time_hours < mod.hrs_in_tmp[tmp]: + tmp_to_check = tmp else: + # Otherwise, we check the following timepoints # First we'll check the next timepoint of the starting timepoint and # start with the duration of the starting timepoint tmp_to_check = mod.next_tmp[tmp, mod.water_system_balancing_type] @@ -137,14 +143,15 @@ def determine_arrival_timepoint(mod, tmp, travel_time_hours): while hours_from_departure_tmp < travel_time_hours: # If we haven't exceeded the travel time yet, we move on to the next tmp # In a 'linear' horizon setting, once we reach the last - # timepoint of the horizon, we break out of the loop since there - # are no more timepoints to consider + # timepoint of the horizon, we set the arrival timepoint to + # "tmp_outside_horizon" and break out of the loop if check_if_boundary_type_and_last_timepoint( mod=mod, tmp=tmp_to_check, balancing_type=mod.water_system_balancing_type, boundary_type="linear", ): + tmp_to_check = "tmp_outside_horizon" break # In a 'circular' horizon setting, once we reach timepoint *t*, # we break out of the loop since there are no more timepoints to diff --git a/gridpath/system/water/water_node_balance.py b/gridpath/system/water/water_node_balance.py index 808ee0e15..4c7add728 100644 --- a/gridpath/system/water/water_node_balance.py +++ b/gridpath/system/water/water_node_balance.py @@ -59,7 +59,7 @@ def add_model_components( """ # ### Expressions ### # - def gross_node_inflow_rate_init(mod, wn, tmp): + def node_inflow_rate_init(mod, wn, tmp): """ Exogenous inflow to node + sum of flow on all links to note in the timepoint of arrival @@ -73,30 +73,24 @@ def gross_node_inflow_rate_init(mod, wn, tmp): m.Gross_Water_Node_Inflow_Rate_Vol_Per_Sec = Expression( m.WATER_NODES, m.TMPS, - initialize=gross_node_inflow_rate_init, + initialize=node_inflow_rate_init, ) - def gross_node_release_rate_vol_per_sec(mod, wn, tmp): - # If we have e reservoir, this is controlled by the reservoir - # variables; otherwise, just set it to inflow for the mass balance - return ( - mod.Gross_Reservoir_Release_Rate_Vol_Per_Sec[wn, tmp] - if wn in mod.WATER_NODES_W_RESERVOIRS - else 0 - ) + ( - mod.Gross_Water_Node_Inflow_Rate_Vol_Per_Sec[wn, tmp] - if wn not in mod.WATER_NODES_W_RESERVOIRS - else 0 + def node_outflow_rate_init(mod, wn, tmp): + return sum( + mod.Water_Link_Flow_Rate_Vol_per_Sec[wl, dep_tmp, arr_tmp] + for (wl, dep_tmp, arr_tmp) in mod.WATER_LINK_DEPARTURE_ARRIVAL_TMPS + if wl in mod.WATER_LINKS_FROM_BY_WATER_NODE[wn] and dep_tmp == tmp ) - m.Gross_Water_Node_Release_Rate_Vol_per_Sec = Expression( + m.Gross_Water_Node_Outflow_Rate_Vol_per_Sec = Expression( m.WATER_NODES, m.TMPS, - initialize=gross_node_release_rate_vol_per_sec, + initialize=node_outflow_rate_init, ) # ### Constraints ### # - def get_total_inflow_volunit(mod, wn, tmp): + def get_total_inflow_for_reservoir_tracking_volunit(mod, wn, tmp): """ Total inflow is exogenous inflow at node plus sum of endogenous inflow from all links to node @@ -109,106 +103,110 @@ def get_total_inflow_volunit(mod, wn, tmp): return inflow_in_tmp - def get_total_release_volunit(mod, wn, tmp): + def get_total_reservoir_release_volunit(mod, wn, tmp): outflow_in_tmp = ( - mod.Gross_Water_Node_Release_Rate_Vol_per_Sec[wn, tmp] + mod.Gross_Reservoir_Release_Rate_Vol_Per_Sec[wn, tmp] * 3600 * mod.hrs_in_tmp[tmp] ) return outflow_in_tmp - def water_mass_balance_constraint_rule(mod, wn, tmp): + def reservoir_storage_tracking_rule(mod, wn, tmp): """ """ - # If no reservoir, simply set total inflow equal to total release - if wn not in mod.WATER_NODES_W_RESERVOIRS: - return get_total_inflow_volunit(mod, wn, tmp) == get_total_release_volunit( - mod, wn, tmp + # No constraint in the first timepoint of a linear horizon (no + # previous timepoints for tracking reservoir levels) + if check_if_first_timepoint( + mod=mod, tmp=tmp, balancing_type=mod.water_system_balancing_type + ) and check_boundary_type( + mod=mod, + tmp=tmp, + balancing_type=mod.water_system_balancing_type, + boundary_type="linear", + ): + return Constraint.Skip + # TODO: add linked horizons + elif check_if_first_timepoint( + mod=mod, tmp=tmp, balancing_type=mod.water_system_balancing_type + ) and check_boundary_type( + mod=mod, + tmp=tmp, + balancing_type=mod.water_system_balancing_type, + boundary_type="linked", + ): + current_tmp_starting_water_volume = None + prev_tmp_starting_water_volume = None + + # Inflows and releases + prev_tmp_inflow = None + prev_tmp_outflow = None + raise ( + UserWarning( + "Linked horizons have not been implemented for " + "water system feature." + ) ) - # If the node does have a reservoir, we'll track the water in storage else: - # No constraint in the first timepoint of a linear horizon (no - # previous timepoint for tracking) - if check_if_first_timepoint( - mod=mod, tmp=tmp, balancing_type=mod.water_system_balancing_type - ) and check_boundary_type( - mod=mod, - tmp=tmp, - balancing_type=mod.water_system_balancing_type, - boundary_type="linear", - ): - return Constraint.Skip - else: - # TODO: add linked horizons - if check_if_first_timepoint( - mod=mod, tmp=tmp, balancing_type=mod.water_system_balancing_type - ) and check_boundary_type( - mod=mod, - tmp=tmp, - balancing_type=mod.water_system_balancing_type, - boundary_type="linked", - ): - current_tmp_starting_water_volume = None - prev_tmp_starting_water_volume = None - - # Inflows and releases - prev_tmp_inflow = None - prev_tmp_release = None - raise ( - UserWarning( - "Linked horizons have not been implemented for " - "water system feature." - ) - ) - else: - current_tmp_starting_water_volume = ( - mod.Reservoir_Starting_Volume_WaterVolumeUnit[wn, tmp] - ) - prev_tmp_starting_water_volume = ( - mod.Reservoir_Starting_Volume_WaterVolumeUnit[ - wn, mod.prev_tmp[tmp, mod.water_system_balancing_type] - ] - ) - - # Inflows and releases; these are already calculated - # based on per sec flows and hours in the timepoint - prev_tmp_inflow = get_total_inflow_volunit( - mod, wn, mod.prev_tmp[tmp, mod.water_system_balancing_type] - ) - - prev_tmp_release = get_total_release_volunit( - mod, wn, mod.prev_tmp[tmp, mod.water_system_balancing_type] - ) - - return current_tmp_starting_water_volume == ( - prev_tmp_starting_water_volume + prev_tmp_inflow - prev_tmp_release - ) + current_tmp_starting_water_volume = ( + mod.Reservoir_Starting_Volume_WaterVolumeUnit[wn, tmp] + ) + prev_tmp_starting_water_volume = ( + mod.Reservoir_Starting_Volume_WaterVolumeUnit[ + wn, mod.prev_tmp[tmp, mod.water_system_balancing_type] + ] + ) - m.Water_Mass_Balance_Constraint = Constraint( - m.WATER_NODES, m.TMPS, rule=water_mass_balance_constraint_rule + # Inflows and releases; these are already calculated + # based on per sec flows and hours in the timepoint + prev_tmp_inflow = get_total_inflow_for_reservoir_tracking_volunit( + mod, wn, mod.prev_tmp[tmp, mod.water_system_balancing_type] + ) + + prev_tmp_outflow = get_total_reservoir_release_volunit( + mod, wn, mod.prev_tmp[tmp, mod.water_system_balancing_type] + ) + + return current_tmp_starting_water_volume == ( + prev_tmp_starting_water_volume + prev_tmp_inflow - prev_tmp_outflow + ) + + m.Reservoir_Storage_Tracking_Constraint = Constraint( + m.WATER_NODES_W_RESERVOIRS, m.TMPS, rule=reservoir_storage_tracking_rule ) # Set the sum of outflows from the node to be equal to discharge & spill - def enforce_outflow_rule(mod, wn, tmp): + # for nodes with storage and to inflow for nodes without storage + def enforce_mass_balance_outflow_rule(mod, wn, tmp): """ - The sum of the flows on all links from this node must equal the gross - outflow from the node. Skip constraint for the last node in the + The sum of the flows on all links from this node must equal the + reservoir release for nodes with reservoirs and total inflows for + reservoirs without reservoirs. Skip constraint for the last node in the network with no out links. + + For linear horizons, the lwater outflows may arrive outside of the + horizon boundary if travel time is more than hours in the remaining + timepoints. We still need to enforce outflow constraints (that are + based on the departure timepoint). These flows have the + "tmp_outside_horizon" index for the arrival timepoint. """ if [wl for wl in mod.WATER_LINKS_FROM_BY_WATER_NODE[wn]]: - return ( - sum( - mod.Water_Link_Flow_Rate_Vol_per_Sec[wl, dep_tmp, arr_tmp] - for (wl, dep_tmp, arr_tmp) in mod.WATER_LINK_DEPARTURE_ARRIVAL_TMPS - if wl in mod.WATER_LINKS_FROM_BY_WATER_NODE[wn] and dep_tmp == tmp + # For nodes with reservoirs, set to reservoir release + if wn in mod.WATER_NODES_W_RESERVOIRS: + return ( + mod.Gross_Water_Node_Outflow_Rate_Vol_per_Sec[wn, tmp] + == mod.Gross_Reservoir_Release_Rate_Vol_Per_Sec[wn, tmp] + ) + else: + # For nodes without reservoirs, set to inflow + return ( + mod.Gross_Water_Node_Outflow_Rate_Vol_per_Sec[wn, tmp] + == mod.Gross_Water_Node_Inflow_Rate_Vol_Per_Sec[wn, tmp] ) - == mod.Gross_Water_Node_Release_Rate_Vol_per_Sec[wn, tmp] - ) else: return Constraint.Skip m.Water_Node_Outflow_Constraint = Constraint( - m.WATER_NODES, m.TMPS, rule=enforce_outflow_rule + m.WATER_NODES, m.TMPS, rule=enforce_mass_balance_outflow_rule ) @@ -328,7 +326,7 @@ def export_results( if wn in m.WATER_NODES_W_RESERVOIRS else None ), - value(m.Gross_Water_Node_Release_Rate_Vol_per_Sec[wn, tmp]), + value(m.Gross_Water_Node_Outflow_Rate_Vol_per_Sec[wn, tmp]), ] for wn in m.WATER_NODES for tmp in m.TMPS