Optionally allow slack flows in first timepoint of linear horizons (#1196)

This can help avoid infeasibilities in the first timepoint of linear horizons
 when flows have not yet reached a particular node, making it impossible to meet
 a minimum flow constraint. This problem could arise in later timepoints too;
 slack is currently only allowed in the first timepoint of the horizon (if enabled
 by the user).

Author: anamileva

db/csvs_test_examples/water/water_system_params/1_base_water_system_params.csv

@@ -1,2 +1,2 @@
-water_system_balancing_type,theoretical_power_coefficient
-day,84.7
+water_system_balancing_type,theoretical_power_coefficient,allow_lin_hrz_first_tmp_flow_slack,allow_lin_hrz_first_tmp_flow_slack_tuning_cost
+day,84.7,0,0

db/db_schema.sql

@@ -1057,9 +1057,11 @@ CREATE TABLE subscenarios_system_water_system_params
 DROP TABLE IF EXISTS inputs_system_water_system_params;
 CREATE TABLE inputs_system_water_system_params
 (
-    water_system_params_scenario_id INTEGER PRIMARY KEY,
-    water_system_balancing_type     TEXT,
-    theoretical_power_coefficient   FLOAT,
+    water_system_params_scenario_id                INTEGER PRIMARY KEY,
+    water_system_balancing_type                    TEXT,
+    theoretical_power_coefficient                  FLOAT,
+    allow_lin_hrz_first_tmp_flow_slack             FLOAT,
+    allow_lin_hrz_first_tmp_flow_slack_tuning_cost FLOAT,
     FOREIGN KEY (water_system_params_scenario_id) REFERENCES
         subscenarios_system_water_system_params (water_system_params_scenario_id)
 );
@@ -6856,16 +6858,17 @@ DROP TABLE IF EXISTS results_system_ra;
 DROP TABLE IF EXISTS results_system_water_link_timepoint;
 CREATE TABLE results_system_water_link_timepoint
 (
-    scenario_id            INTEGER,
-    weather_iteration      INTEGER,
-    hydro_iteration        INTEGER,
-    availability_iteration INTEGER,
-    subproblem_id          INTEGER,
-    stage_id               INTEGER,
-    water_link             VARCHAR(32),
-    departure_timepoint    INTEGER,
-    arrival_timepoint      INTEGER,
-    water_flow_vol_per_sec FLOAT,
+    scenario_id                       INTEGER,
+    weather_iteration                 INTEGER,
+    hydro_iteration                   INTEGER,
+    availability_iteration            INTEGER,
+    subproblem_id                     INTEGER,
+    stage_id                          INTEGER,
+    water_link                        VARCHAR(32),
+    departure_timepoint               INTEGER,
+    arrival_timepoint                 INTEGER,
+    water_flow_vol_per_sec            FLOAT,
+    water_flow_slack_used_vol_per_sec FLOAT,
     PRIMARY KEY (scenario_id, weather_iteration, hydro_iteration,
                  availability_iteration, subproblem_id, stage_id, water_link,
                  departure_timepoint)
@@ -6973,6 +6976,7 @@ CREATE TABLE results_system_costs
     Total_Carbon_Credit_Costs                               FLOAT,
     Total_Peak_Deviation_Monthly_Demand_Charge_Cost         FLOAT,
     Total_Policy_Target_Balance_Penalty_Costs               FLOAT,
+    Total_Flow_Constraint_Slack_Tuning_Cost                 FLOAT,
     PRIMARY KEY (scenario_id, weather_iteration, hydro_iteration,
                  availability_iteration, subproblem_id, stage_id)
 );

examples/hydro_system_exog_elev/hydro_iteration_1/inputs/water_system_params.tab

@@ -1,2 +1,2 @@
-water_system_balancing_type	theoretical_power_coefficient
-day	84.7
+water_system_balancing_type	theoretical_power_coefficient	allow_lin_hrz_first_tmp_flow_slack	allow_lin_hrz_first_tmp_flow_slack_tuning_cost
+day	84.7	0.0	0.0

examples/hydro_system_exog_elev_w_travel_time/hydro_iteration_1/inputs/water_system_params.tab

@@ -1,2 +1,2 @@
-water_system_balancing_type	theoretical_power_coefficient
-day	84.7
+water_system_balancing_type	theoretical_power_coefficient	allow_lin_hrz_first_tmp_flow_slack	allow_lin_hrz_first_tmp_flow_slack_tuning_cost
+day	84.7	0.0	0.0

gridpath/auxiliary/module_list.py

@@ -89,8 +89,8 @@ def all_modules_list():
         "system.reliability.local_capacity.local_capacity_requirement",
         "system.markets.prices",
         "system.water.water_system_params",
-        "system.water.water_flows",
         "system.water.water_nodes",
+        "system.water.water_flows",
         "system.water.reservoirs",
         "system.water.water_node_balance",
         "system.water.powerhouses",
@@ -266,6 +266,7 @@ def all_modules_list():
         "objective.system.reliability.local_capacity"
         ".aggregate_local_capacity_violation_penalties",
         "objective.system.aggregate_market_revenue_and_costs",
+        "objective.system.water.aggregate_flow_constraint_slack_use_tuning_costs",
         "objective.max_npv",
     ]
     return all_modules
@@ -491,6 +492,7 @@ def optional_modules_list():
             "system.water.water_node_balance",
             "system.water.water_flows",
             "system.water.powerhouses",
+            "objective.system.water.aggregate_flow_constraint_slack_use_tuning_costs",
         ],
         "tuning": [
             "project.operations.tuning_costs",

gridpath/objective/system/water/__init__.py

@@ -0,0 +1,71 @@
+# Copyright 2016-2024 Blue Marble Analytics LLC.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+
+from pyomo.environ import Param, Expression
+
+from gridpath.auxiliary.dynamic_components import cost_components
+
+
+def add_model_components(
+    m,
+    d,
+    scenario_directory,
+    weather_iteration,
+    hydro_iteration,
+    availability_iteration,
+    subproblem,
+    stage,
+):
+    """
+
+    :param m:
+    :param d:
+    :return:
+    """
+
+    def total_flow_constraint_slack_tuning_cost_rule(mod):
+        """
+        Ramp tuning costs for all projects
+        :param mod:
+        :return:
+        """
+        return sum(
+            mod.Water_Link_Flow_Rate_Vol_per_Sec_LinHrzFirstTmpSlack[wl, tmp]
+            * mod.allow_lin_hrz_first_tmp_flow_slack_tuning_cost
+            * mod.hrs_in_tmp[tmp]
+            * mod.tmp_weight[tmp]
+            * mod.number_years_represented[mod.period[tmp]]
+            * mod.discount_factor[mod.period[tmp]]
+            for wl in mod.WATER_LINKS
+            for tmp in mod.TMPS
+        )
+
+    m.Total_Flow_Constraint_Slack_Tuning_Cost = Expression(
+        rule=total_flow_constraint_slack_tuning_cost_rule
+    )
+
+    record_dynamic_components(dynamic_components=d)
+
+
+def record_dynamic_components(dynamic_components):
+    """
+    :param dynamic_components:
+
+    Add tuning costs to cost components
+    """
+
+    getattr(dynamic_components, cost_components).append(
+        "Total_Flow_Constraint_Slack_Tuning_Cost"
+    )

gridpath/objective/system/water/aggregate_flow_constraint_slack_use_tuning_costs.py

@@ -34,6 +34,7 @@
 from gridpath.project.common_functions import (
     check_if_boundary_type_and_last_timepoint,
     check_boundary_type,
+    check_if_boundary_type_and_first_timepoint,
 )
 
 
@@ -69,7 +70,7 @@ def water_link_departure_arrival_tmp_init(mod):
             for departure_tmp in mod.TMPS:
                 arrival_tmp = determine_arrival_timepoint(
                     mod=mod,
-                    tmp=departure_tmp,
+                    dep_tmp=departure_tmp,
                     travel_time_hours=mod.water_link_flow_transport_time_hours[wl],
                 )
                 if arrival_tmp is not None:
@@ -89,17 +90,48 @@ def water_link_departure_arrival_tmp_init(mod):
         m.WATER_LINK_DEPARTURE_ARRIVAL_TMPS, within=NonNegativeReals
     )
 
+    m.Water_Link_Flow_Rate_Vol_per_Sec_LinHrzFirstTmpSlack = Var(
+        m.WATER_LINKS, m.TMPS, within=NonNegativeReals, initialize=0
+    )
+
     # ### Constraints ### #
     def min_flow_rule(mod, wl, dep_tmp, arr_tmp):
         return (
             mod.Water_Link_Flow_Rate_Vol_per_Sec[wl, dep_tmp, arr_tmp]
+            + mod.Water_Link_Flow_Rate_Vol_per_Sec_LinHrzFirstTmpSlack[wl, dep_tmp]
             >= mod.min_flow_vol_per_second[wl, dep_tmp]
         )
 
     m.Water_Link_Minimum_Flow_Constraint = Constraint(
         m.WATER_LINK_DEPARTURE_ARRIVAL_TMPS, rule=min_flow_rule
     )
 
+    def min_flow_slack_rule(mod, wl, tmp):
+        if mod.allow_lin_hrz_first_tmp_flow_slack:
+            if check_if_boundary_type_and_first_timepoint(
+                mod=mod,
+                tmp=tmp,
+                boundary_type="linear",
+                balancing_type=mod.water_system_balancing_type,
+            ):
+                return (
+                    mod.Water_Link_Flow_Rate_Vol_per_Sec_LinHrzFirstTmpSlack[wl, tmp]
+                    <= mod.min_flow_vol_per_second[wl, tmp]
+                )
+            else:
+                return (
+                    mod.Water_Link_Flow_Rate_Vol_per_Sec_LinHrzFirstTmpSlack[wl, tmp]
+                    == 0
+                )
+        else:
+            return (
+                mod.Water_Link_Flow_Rate_Vol_per_Sec_LinHrzFirstTmpSlack[wl, tmp] == 0
+            )
+
+    m.MinFlowSlackConstraint = Constraint(
+        m.WATER_LINKS, m.TMPS, rule=min_flow_slack_rule
+    )
+
     def max_flow_rule(mod, wl, dep_tmp, arr_tmp):
         return (
             mod.Water_Link_Flow_Rate_Vol_per_Sec[wl, dep_tmp, arr_tmp]
@@ -111,60 +143,61 @@ def max_flow_rule(mod, wl, dep_tmp, arr_tmp):
     )
 
 
-def determine_arrival_timepoint(mod, tmp, travel_time_hours):
+def determine_arrival_timepoint(mod, dep_tmp, travel_time_hours):
+    # If travel time is less than the hours in the departure timepoint, we stay
+    # in the same timepoint
+    if travel_time_hours < mod.hrs_in_tmp[dep_tmp]:
+        arr_tmp = dep_tmp
     # 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(
+    # timepoints to check and we'll return 'tmp_outside_horizon'
+    elif check_if_boundary_type_and_last_timepoint(
         mod=mod,
-        tmp=tmp,
+        tmp=dep_tmp,
         balancing_type=mod.water_system_balancing_type,
         boundary_type="linear",
     ):
-        tmp_to_check = "tmp_outside_horizon"
+        arr_tmp = "tmp_outside_horizon"
     elif check_if_boundary_type_and_last_timepoint(
         mod=mod,
-        tmp=tmp,
+        tmp=dep_tmp,
         balancing_type=mod.water_system_balancing_type,
         boundary_type="linked",
     ):
         # TODO: add linked
-        tmp_to_check = None
-    # 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
+        arr_tmp = None
     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]
-        hours_from_departure_tmp = mod.hrs_in_tmp[tmp]
+        arr_tmp = mod.next_tmp[dep_tmp, mod.water_system_balancing_type]
+        hours_from_departure_tmp = mod.hrs_in_tmp[dep_tmp]
         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 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,
+                tmp=arr_tmp,
                 balancing_type=mod.water_system_balancing_type,
                 boundary_type="linear",
             ):
-                tmp_to_check = "tmp_outside_horizon"
+                arr_tmp = "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
-            # consider (we have already checked all horizon timepoints)
+            # In a 'circular' horizon setting, once we loop back to the
+            # departure timepoint again, we break out of the loop since there
+            # are no more timepoints to consider (we have already checked all
+            # horizon timepoints)
             elif (
                 check_boundary_type(
                     mod=mod,
-                    tmp=tmp,
+                    tmp=dep_tmp,
                     balancing_type=mod.water_system_balancing_type,
                     boundary_type="circular",
                 )
-                and tmp_to_check == tmp
+                and arr_tmp == dep_tmp
             ):
+                arr_tmp = "tmp_outside_horizon"
                 break
             # TODO: only allow the first horizon of a subproblem to have
             #  linked timepoints
@@ -173,21 +206,20 @@ def determine_arrival_timepoint(mod, tmp, travel_time_hours):
             # timepoints until we reach the target min time
             elif check_if_boundary_type_and_last_timepoint(
                 mod=mod,
-                tmp=tmp_to_check,
+                tmp=arr_tmp,
                 balancing_type=mod.water_system_balancing_type,
                 boundary_type="linked",
             ):
                 # TODO: add linked
+                arr_tmp = None
                 break
             # Otherwise, we move on to the next timepoint and will add that
             # timepoint's duration to hours_from_departure_tmp
             else:
-                hours_from_departure_tmp += mod.hrs_in_tmp[tmp_to_check]
-                tmp_to_check = mod.next_tmp[
-                    tmp_to_check, mod.water_system_balancing_type
-                ]
+                hours_from_departure_tmp += mod.hrs_in_tmp[arr_tmp]
+                arr_tmp = mod.next_tmp[arr_tmp, mod.water_system_balancing_type]
 
-    return tmp_to_check
+    return arr_tmp
 
 
 def load_model_data(
@@ -384,13 +416,15 @@ def export_results(
     """
     results_columns = [
         "water_flow_vol_per_sec",
+        "water_flow_slack_used_vol_per_sec",
     ]
     data = [
         [
             wl,
             dep_tmp,
             arr_tmp,
             value(m.Water_Link_Flow_Rate_Vol_per_Sec[wl, dep_tmp, arr_tmp]),
+            value(m.Water_Link_Flow_Rate_Vol_per_Sec_LinHrzFirstTmpSlack[wl, dep_tmp]),
         ]
         for (wl, dep_tmp, arr_tmp) in m.WATER_LINK_DEPARTURE_ARRIVAL_TMPS
     ]

gridpath/system/water/water_flows.py

@@ -210,33 +210,6 @@ def enforce_mass_balance_outflow_rule(mod, wn, tmp):
     )
 
 
-def load_model_data(
-    m,
-    d,
-    data_portal,
-    scenario_directory,
-    weather_iteration,
-    hydro_iteration,
-    availability_iteration,
-    subproblem,
-    stage,
-):
-
-    data_portal.load(
-        filename=os.path.join(
-            scenario_directory,
-            weather_iteration,
-            hydro_iteration,
-            availability_iteration,
-            subproblem,
-            stage,
-            "inputs",
-            "water_inflows.tab",
-        ),
-        param=m.exogenous_water_inflow_rate_vol_per_sec,
-    )
-
-
 def validate_inputs(
     scenario_id,
     subscenarios,