diff --git a/gridpath/auxiliary/scaling.py b/gridpath/auxiliary/scaling.py new file mode 100644 index 0000000000..41e72eeead --- /dev/null +++ b/gridpath/auxiliary/scaling.py @@ -0,0 +1,366 @@ +# 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. + +""" +Numerical scaling for improved solver conditioning. + +Large models can span a wide numerical range: power/energy quantities in MW/MWh +(``10**0``-``10**5``) sit in the same matrix as dollar penalties in ``$/MWh`` +that reach ``10**6`` or more after net-present-value weighting. Solvers warn +about (and can stall on) such coefficient ranges. + +This module assigns a Pyomo ``scaling_factor`` Suffix so that the built model +can be solved in scaled units -- e.g. GW/GWh instead of MW/MWh (a power/energy +scale factor of ``1000``) and millions of dollars instead of dollars (a dollar +scale factor of ``1000000``). The actual reformulation and the inverse mapping +of the solution (variable values, duals, and reduced costs) back to native +units are done by Pyomo's ``core.scale_model`` transformation; this module only +decides the per-component factors. + +Why this is safe: ``core.scale_model`` is an *exact* affine reformulation for +any positive factor assignment (a component with no factor defaults to ``1.0``, +i.e. identity). A units misclassification here therefore only affects +conditioning -- it can never change the optimal solution. In particular, the +mixed-unit cost coefficients (``$/MWh``, ``$/MW-yr``, ``$/MMBtu``, ...) live as +*constants inside rows*, not as variables, and the transformation rescales every +constant coefficient in every row automatically; they are handled structurally +and never need to be enumerated here. + +The classification of a component's units is therefore a conditioning heuristic +driven by GridPath's naming conventions: + + * variables carrying an ``MW`` / ``MWh`` / ``MWs`` token (or a trailing + ``Power``) are power/energy; + * variables carrying a ``Cost`` token are dollars; + * integer/binary variables are never scaled (scaling their bounds would + break integrality); + * a constraint's factor is inferred from the variables in its body; + * the objective (net present value, in dollars) is scaled by the dollar + factor. + +Anything not recognized is left at ``1.0`` (unscaled) -- correct, if not +maximally conditioned. Separate "commodity" chains (fuel in MMBtu, emissions in +tons, water volumes) fall into this bucket; extending the heuristic to those is +a localized change here if a model needs it. +""" + +from pyomo.environ import Suffix, Var, Constraint, Objective +from pyomo.core.expr import identify_variables, replace_expressions + + +# Name tokens (split on "_") that mark a variable as power/energy or dollars. +POWER_ENERGY_TOKENS = frozenset({"MW", "MWh", "MWs"}) +DOLLAR_TOKEN = "Cost" + + +def classify_variable_units(name): + """Classify a variable by its name using GridPath naming conventions. + + Args: + name: The variable component's name (e.g. ``"GenSimple_Provide_Power_MW"``). + + Returns: + One of ``"power"``, ``"dollar"``, or ``None`` (unrecognized -- leave + unscaled). Dollars take precedence over power, though no GridPath + variable name currently carries both tokens. + """ + tokens = name.split("_") + if DOLLAR_TOKEN in tokens: + return "dollar" + if any(t in POWER_ENERGY_TOKENS for t in tokens): + return "power" + # A trailing "Power" token catches unsuffixed power variables such as + # "Net_Market_Purchased_Power" without matching generic-unit names like + # "Fuel_Prod_Consume_Power_PowerUnit" (whose last token is "PowerUnit"). + if tokens[-1] == "Power": + return "power" + return None + + +def _variable_factor(var, s_power, s_dollar): + """Return the scaling factor for a variable container. + + Integer/binary variables are never scaled (factor ``1.0``): the scaling + transformation only rescales bounds and values, not the domain, so scaling + an integer variable by a non-integer factor would silently break + integrality. + + Args: + var: A Pyomo ``Var`` container. + s_power: The power/energy scaling factor (``1 / power_scale_factor``). + s_dollar: The dollar scaling factor (``1 / dollar_scale_factor``). + + Returns: + The float scaling factor to assign to the variable (``1.0`` if it should + not be scaled). + """ + representative = next(iter(var.values()), None) + if representative is None: + return 1.0 + if representative.is_integer(): + return 1.0 + + units = classify_variable_units(var.name) + if units == "power": + return s_power + if units == "dollar": + return s_dollar + return 1.0 + + +def _constraint_factor(constraint, container_factor, s_dollar): + """Infer a constraint's scaling factor from the variables in its body. + + A GridPath ``Constraint`` container is built from a single rule over an + index set, so every data object shares the same variable-unit structure; we + therefore inspect only one representative active data object rather than the + whole (potentially huge) index set. + + The rule, applied to the set of non-unity variable factors appearing in the + body: + + * none -> ``1.0`` (nothing to condition); + * the dollar factor appears -> use the dollar factor. This is a + cost-definition row (e.g. ``Hurdle_Cost >= flow * rate``); scaling by + the dollar factor makes the defining dollar variable's coefficient + exactly ``1``; + * exactly one factor -> use it. This is a homogeneous row (e.g. a power + balance ``sum(MW) == load``); coefficients stay ``O(1)`` and the + constant right-hand side is rescaled into the new unit; + * more than one non-dollar factor -> use the smallest (largest scale), + a safe fallback for the rare genuinely mixed row. + + Args: + constraint: A Pyomo ``Constraint`` container. + container_factor: Map from ``id`` of a ``Var`` container to its assigned + factor. + s_dollar: The dollar scaling factor. + + Returns: + The float scaling factor to assign to the constraint. + """ + representative = next((cd for cd in constraint.values() if cd.active), None) + if representative is None: + return 1.0 + + factors = set() + for var_data in identify_variables(representative.body, include_fixed=False): + factor = container_factor.get(id(var_data.parent_component()), 1.0) + factors.add(factor) + factors.discard(1.0) + + if not factors: + return 1.0 + if s_dollar in factors: + return s_dollar + return min(factors) + + +def assign_scaling_factors(instance, power_scale_factor, dollar_scale_factor): + """Attach a ``scaling_factor`` Suffix to a built model instance. + + The suffix is consumed by Pyomo's ``core.scale_model`` transformation. Power + and energy quantities are divided by ``power_scale_factor`` (e.g. ``1000`` + for MW->GW, MWh->GWh) and dollar quantities by ``dollar_scale_factor`` (e.g. + ``1000000`` for $->$M). Factors are assigned at the component-container level + (one suffix entry per container, not per data object), which the + transformation's suffix lookup resolves to every contained data object -- + essential to keep the suffix small on large models. + + Args: + instance: A concrete Pyomo model instance (already built and, if + applicable, with fixed variables set). + power_scale_factor: Factor to divide power/energy quantities by. Must be + positive. + dollar_scale_factor: Factor to divide dollar quantities by. Must be + positive. + + Returns: + The instance, with an ``instance.scaling_factor`` Suffix attached. + + Raises: + ValueError: If either scale factor is not positive. + """ + if power_scale_factor <= 0 or dollar_scale_factor <= 0: + raise ValueError( + "Scale factors must be positive; got power_scale_factor=" + f"{power_scale_factor}, dollar_scale_factor={dollar_scale_factor}." + ) + + s_power = 1.0 / power_scale_factor + s_dollar = 1.0 / dollar_scale_factor + + instance.scaling_factor = Suffix(direction=Suffix.EXPORT) + + # Variables: classify each container and record its factor for the + # constraint inference pass. + container_factor = {} + for var in instance.component_objects(Var, descend_into=True): + factor = _variable_factor(var, s_power, s_dollar) + container_factor[id(var)] = factor + if factor != 1.0: + instance.scaling_factor[var] = factor + + # Constraints: infer each container's factor from the variables in its body. + for constraint in instance.component_objects( + Constraint, descend_into=True, active=True + ): + factor = _constraint_factor(constraint, container_factor, s_dollar) + if factor != 1.0: + instance.scaling_factor[constraint] = factor + + # Objective(s): net present value is in dollars. + for objective in instance.component_objects( + Objective, descend_into=True, active=True + ): + instance.scaling_factor[objective] = s_dollar + + return instance + + +def propagate_scaled_solution(scaled_instance, instance): + """Map a solved scaled model's solution back onto the original instance. + + This mirrors Pyomo's ``ScaleModel.propagate_solution`` (variable values are + divided by their scaling factor; duals are multiplied by the constraint + factor and divided by the objective factor) but tolerates constraints for + which the solver did not return a dual. Pyomo's own implementation assumes + every constraint has a dual and raises ``KeyError`` otherwise; solvers + routinely leave some duals unpopulated (e.g. non-binding market limits under + CBC), and GridPath's export path already handles a missing dual as ``None``, + so skipping them here keeps the scaled and unscaled paths equivalent. + + The ``scaled_instance`` must have been produced by the ``core.scale_model`` + transformation (it carries the ``component_scaling_factor_map`` and + ``scaled_component_to_original_name_map`` used for back-mapping). + + Args: + scaled_instance: The scaled model, after it has been solved. + instance: The original (native-unit) model to receive the solution. + + Returns: + The original ``instance`` with variable values and duals populated in + native units. + """ + factor_map = scaled_instance.component_scaling_factor_map + name_map = scaled_instance.scaled_component_to_original_name_map + + has_dual = hasattr(scaled_instance, "dual") and hasattr(instance, "dual") + + # Objective scaling factor (duals/reduced costs are relative to it). There + # is exactly one active objective in GridPath (the NPV). + objective_factor = 1.0 + for scaled_obj in scaled_instance.component_data_objects( + Objective, active=True, descend_into=True + ): + objective_factor = factor_map[scaled_obj] + break + + # Variable values: original = scaled / factor. + for scaled_var in scaled_instance.component_objects(Var, descend_into=True): + original_var = instance.find_component(name_map[scaled_var]) + for k in scaled_var: + scaled_value = scaled_var[k].value + if scaled_value is None: + original_var[k].set_value(None, skip_validation=True) + else: + original_var[k].set_value( + scaled_value / factor_map[scaled_var[k]], + skip_validation=True, + ) + + # Duals: original = scaled * constraint_factor / objective_factor. Skip any + # constraint the solver left without a dual. + if has_dual: + for scaled_con in scaled_instance.component_objects( + Constraint, descend_into=True + ): + original_con = instance.find_component(name_map[scaled_con]) + for k in scaled_con: + if scaled_con[k] not in scaled_instance.dual: + continue + instance.dual[original_con[k]] = ( + scaled_instance.dual[scaled_con[k]] + * factor_map[scaled_con[k]] + / objective_factor + ) + + return instance + + +def invert_scaled_solution_in_place(instance): + """Un-scale a solved model that was scaled in place, restoring native units. + + Counterpart to ``propagate_scaled_solution`` for the in-place path: instead + of cloning the model (``create_using``) and mapping the solution back onto a + pristine original, the model itself was scaled with ``apply_to(rename=False)`` + and solved. This reverses the scaling on the solved model so that everything + downstream reads native units: + + * Objective expression: ``apply_to`` rewrote it as + ``s_obj * (expression in scaled variables)``. We substitute each scaled + variable ``v -> v * s_v`` (recovering the native-variable expression, + still multiplied by ``s_obj``) and then divide the whole objective by + ``s_obj``. This matters because ``save_objective_function_value`` reads + ``instance.NPV()`` directly; without this it would be off by ``s_obj``. + * Variable values: ``v <- v / s_v`` (native units). + * Duals: ``dual <- dual * s_c / s_obj`` (same formula as + ``propagate_scaled_solution``). + + Named cost/revenue ``Expression`` components are NOT rewritten by scaling + (only ``Constraint``/``Objective``/``Var`` are), so once variable values are + native they already evaluate to native dollars and need no adjustment. + + The instance must have been scaled with ``TransformationFactory( + 'core.scale_model').apply_to(instance, rename=False)`` (which stores the + ``component_scaling_factor_map`` this reads) and then solved. + + Args: + instance: The in-place-scaled, solved model to restore to native units. + + Returns: + The instance, with objective, variable values, and duals in native units. + """ + factor_map = instance.component_scaling_factor_map + + objective_factor = 1.0 + active_objectives = list( + instance.component_data_objects(Objective, active=True, descend_into=True) + ) + for obj in active_objectives: + objective_factor = factor_map[obj] + break + + # Reverse the objective substitution: undo v -> v / s_v (i.e. put back + # v -> v * s_v), then remove the objective row factor s_obj. + for obj in active_objectives: + substitution = {id(v): v * factor_map[v] for v in identify_variables(obj.expr)} + obj.expr = replace_expressions(obj.expr, substitution) / objective_factor + + # Un-scale variable values. + for var in instance.component_data_objects(Var, descend_into=True): + if var.value is not None: + var.set_value(var.value / factor_map[var], skip_validation=True) + + # Rescale duals (skip any constraint the solver left without a dual). + if hasattr(instance, "dual"): + for constraint in instance.component_data_objects( + Constraint, active=True, descend_into=True + ): + if constraint in instance.dual: + instance.dual[constraint] = ( + instance.dual[constraint] + * factor_map[constraint] + / objective_factor + ) + + return instance diff --git a/gridpath/common_functions.py b/gridpath/common_functions.py index 7a0ae3895c..8ecfcb779a 100644 --- a/gridpath/common_functions.py +++ b/gridpath/common_functions.py @@ -347,6 +347,35 @@ def get_run_scenario_parser(): help="Skip quick summary text file", ) + # Numerical scaling (for solver conditioning) + parser.add_argument( + "--power_scale_factor", + default=1.0, + type=float, + help="Divide power/energy quantities (MW, MWh) by this factor when " + "solving, then map the solution back to native units. E.g. 1000 solves " + "in GW/GWh. Default 1.0 (no scaling). See gridpath/auxiliary/scaling.py.", + ) + parser.add_argument( + "--dollar_scale_factor", + default=1.0, + type=float, + help="Divide dollar quantities by this factor when solving, then map " + "the solution back to native units. E.g. 1000000 solves in millions of " + "dollars. Default 1.0 (no scaling). See gridpath/auxiliary/scaling.py.", + ) + parser.add_argument( + "--scale_mode", + default="out_of_place", + choices=["out_of_place", "in_place"], + help="How to apply numerical scaling (only relevant if a scale factor " + "is set). 'out_of_place' (default) solves a scaled clone and maps the " + "solution back, keeping the original model pristine. 'in_place' scales " + "the model itself and inverts the solution afterward, avoiding the clone " + "(~half the peak memory and faster setup) at the cost of mutating the " + "model; use it for very large problems where the clone is expensive.", + ) + return parser diff --git a/gridpath/run_scenario.py b/gridpath/run_scenario.py index 4724d0df95..2eb70ec66a 100644 --- a/gridpath/run_scenario.py +++ b/gridpath/run_scenario.py @@ -38,6 +38,7 @@ SolverFactory, SolverStatus, TerminationCondition, + TransformationFactory, ) # from pyomo.util.infeasible import log_infeasible_constraints @@ -67,6 +68,11 @@ ) from gridpath.auxiliary.dynamic_components import DynamicComponents from gridpath.auxiliary.module_list import determine_modules, load_modules +from gridpath.auxiliary.scaling import ( + assign_scaling_factors, + propagate_scaled_solution, + invert_scaled_solution_in_place, +) def start_step(step, quiet): @@ -246,14 +252,73 @@ def create_problem( def solve_problem(parsed_arguments, instance, timing_summary_file_path=None): # Solve + power_scale_factor = getattr(parsed_arguments, "power_scale_factor", 1.0) + dollar_scale_factor = getattr(parsed_arguments, "dollar_scale_factor", 1.0) + + # No scaling requested: solve the instance directly (the default path -- no + # suffix, no clone, no transformation). + if power_scale_factor == 1.0 and dollar_scale_factor == 1.0: + step_start_time = start_step(step="Solving", quiet=parsed_arguments.quiet) + results = solve(instance, parsed_arguments) + report_step_timing( + step="Solving", + step_start_time=step_start_time, + quiet=parsed_arguments.quiet, + timing_summary_file_path=timing_summary_file_path, + ) + return instance, results + + # Scaling requested: assign scaling factors on the instance, then either + # (out_of_place) solve a scaled clone and map the solution back onto the + # pristine original, or (in_place) scale the instance itself and invert the + # solution afterward. Both leave the instance handed downstream in native + # units, so results export / objective / duals code is unchanged. + assign_scaling_factors( + instance, + power_scale_factor=power_scale_factor, + dollar_scale_factor=dollar_scale_factor, + ) + scaler = TransformationFactory("core.scale_model") + scale_mode = getattr(parsed_arguments, "scale_mode", "out_of_place") + + if scale_mode == "in_place": + # Scale the model itself (no clone -> ~half the peak memory and faster + # setup for large models), solve it, then restore native units on the + # same instance. + scaler.apply_to(instance, rename=False) + step_start_time = start_step(step="Solving", quiet=parsed_arguments.quiet) + results = solve(instance, parsed_arguments) + report_step_timing( + step="Solving", + step_start_time=step_start_time, + quiet=parsed_arguments.quiet, + timing_summary_file_path=timing_summary_file_path, + ) + invert_scaled_solution_in_place(instance) + return instance, results + + # out_of_place (default): solve a scaled clone, then map the solution + # (variable values and duals) back onto the original native-unit instance. + scaled_instance = scaler.create_using(instance) step_start_time = start_step(step="Solving", quiet=parsed_arguments.quiet) - results = solve(instance, parsed_arguments) + results = solve(scaled_instance, parsed_arguments) report_step_timing( step="Solving", step_start_time=step_start_time, quiet=parsed_arguments.quiet, timing_summary_file_path=timing_summary_file_path, ) + # Use our own back-mapping rather than scaler.propagate_solution because the + # latter raises if the solver left any constraint without a dual (which + # happens, e.g. non-binding market limits under CBC); ours skips those, as + # GridPath's export path already treats a missing dual as None. + propagate_scaled_solution(scaled_instance, instance) + + # Release the scaled clone promptly (it doubles peak memory for large + # models); matches the garbage-collection discipline elsewhere in this + # module. + del scaled_instance + gc.collect() return instance, results @@ -1716,6 +1781,31 @@ def main(args=None): # Parse arguments parsed_args = parse_arguments(args) + # Numerical scaling is applied at solve time (see solve_problem), so it is + # incompatible with the paths that skip solving and instead load a solution + # from a file (whose values are in unknown units) or only write the problem + # file (which would be written unscaled). Fail fast rather than silently + # mis-handle these. + scaling_requested = ( + parsed_args.power_scale_factor != 1.0 or parsed_args.dollar_scale_factor != 1.0 + ) + if scaling_requested: + incompatible = [ + flag + for flag in ( + "load_cplex_solution", + "load_gurobi_solution", + "load_highs_solution", + "create_lp_problem_file_only", + ) + if getattr(parsed_args, flag, False) + ] + if incompatible: + raise ValueError( + "--power_scale_factor / --dollar_scale_factor cannot be " + "combined with {}.".format(", ".join("--" + f for f in incompatible)) + ) + scenario_directory = determine_scenario_directory( scenario_location=parsed_args.scenario_location, scenario_name=parsed_args.scenario, diff --git a/gridpath/system/load_balance/load_balance.py b/gridpath/system/load_balance/load_balance.py index 4850d3a7af..3d04dbc660 100644 --- a/gridpath/system/load_balance/load_balance.py +++ b/gridpath/system/load_balance/load_balance.py @@ -146,6 +146,11 @@ def meet_load_rule(mod, z, tmp): m.Meet_Load_Constraint = Constraint(m.LOAD_ZONES, m.TMPS, rule=meet_load_rule) def use_limit_constraint_rule(mod, lz): + # No limit specified (defaults to +inf): skip the constraint entirely + # rather than build a row with an infinite (or huge) RHS, which would + # be a free row that only hurts solver scaling. + if mod.unserved_energy_limit_mwh[lz] == float("inf"): + return Constraint.Skip return ( sum( mod.Unserved_Energy_MW_Expression[lz, tmp] @@ -161,6 +166,9 @@ def use_limit_constraint_rule(mod, lz): ) def max_unserved_load_limit_constraint_rule(mod, lz, tmp): + # No limit specified (defaults to +inf): skip (see above). + if mod.max_unserved_load_limit_mw[lz] == float("inf"): + return Constraint.Skip return ( mod.Unserved_Energy_MW_Expression[lz, tmp] <= mod.max_unserved_load_limit_mw[lz] diff --git a/gridpath/transmission/capacity/capacity.py b/gridpath/transmission/capacity/capacity.py index 8bd5fffe35..674d0d791c 100644 --- a/gridpath/transmission/capacity/capacity.py +++ b/gridpath/transmission/capacity/capacity.py @@ -22,6 +22,7 @@ again depend on the line's *capacity_type*. """ +import math import os.path import pandas as pd from pyomo.environ import Set, Expression, value @@ -79,13 +80,28 @@ def add_model_components( | | :code:`TX_OPR_TMPS` | | | | Two-dimensional set of the transmission lines and their operational | - | timepoints, derived from :code:`TX_OPR_PRDS` and the timepoitns in each | + | timepoints, derived from :code:`TX_OPR_PRDS` and the timepoints in each | | period. | +-------------------------------------------------------------------------+ | | :code:`TX_LINES_OPR_IN_TMP` | | | *Defined over*: :code:`TIMEPOINTS` | | | - | Indexed set of transmission lines operatoinal in each timepoint. | + | Indexed set of transmission lines operational in each timepoint. | + +-------------------------------------------------------------------------+ + | | :code:`TX_OPR_PRDS_W_MIN_LIMIT` | + | | + | Subset of :code:`TX_OPR_PRDS` for line-periods that have a lower flow | + | limit. A capacity type may declare a line-period unconstrained (no | + | limit) via :code:`min_limit_is_unconstrained_rule`; capacity types | + | without that method are always constrained (the default). The | + | operational types build their minimum-flow constraints over this | + | subset, so unconstrained line-periods get no such constraint. | + +-------------------------------------------------------------------------+ + | | :code:`TX_OPR_PRDS_W_MAX_LIMIT` | + | | + | Subset of :code:`TX_OPR_PRDS` for line-periods that have an upper flow | + | limit (analogous to :code:`TX_OPR_PRDS_W_MIN_LIMIT`, via | + | :code:`max_limit_is_unconstrained_rule`). | +-------------------------------------------------------------------------+ | @@ -186,6 +202,47 @@ def tx_max_capacity_rule(mod, tx, p): m.Tx_Max_Capacity_MW = Expression(m.TX_OPR_PRDS, rule=tx_max_capacity_rule) + # Sets of line-periods that have a lower / upper flow limit. A capacity + # type may declare a line-period "unconstrained" (no flow limit) by + # defining min_limit_is_unconstrained_rule / max_limit_is_unconstrained_rule + # and returning True; capacity types without those methods are always + # constrained (the default), so no line is ever silently left unbounded. + # The operational types build their min/max flow constraints over these + # subsets, skipping unconstrained line-periods entirely. + def tx_min_limit_is_unconstrained(mod, tx, p): + cap_type = mod.tx_capacity_type[tx] + module = imported_tx_capacity_modules[cap_type] + if hasattr(module, "min_limit_is_unconstrained_rule"): + return module.min_limit_is_unconstrained_rule(mod, tx, p) + return False + + def tx_max_limit_is_unconstrained(mod, tx, p): + cap_type = mod.tx_capacity_type[tx] + module = imported_tx_capacity_modules[cap_type] + if hasattr(module, "max_limit_is_unconstrained_rule"): + return module.max_limit_is_unconstrained_rule(mod, tx, p) + return False + + m.TX_OPR_PRDS_W_MIN_LIMIT = Set( + dimen=2, + within=m.TX_OPR_PRDS, + initialize=lambda mod: [ + (tx, p) + for (tx, p) in mod.TX_OPR_PRDS + if not tx_min_limit_is_unconstrained(mod, tx, p) + ], + ) + + m.TX_OPR_PRDS_W_MAX_LIMIT = Set( + dimen=2, + within=m.TX_OPR_PRDS, + initialize=lambda mod: [ + (tx, p) + for (tx, p) in mod.TX_OPR_PRDS + if not tx_max_limit_is_unconstrained(mod, tx, p) + ], + ) + # Set Rules ############################################################################### @@ -255,12 +312,17 @@ def export_results( "max_mw", ] + # An unconstrained line-period has an infinite capacity; report it as + # NULL rather than the literal "inf" so the results stay numeric. + def _finite_or_none(v): + return None if math.isinf(v) else v + data = [ [ tx_line, prd, - value(m.Tx_Min_Capacity_MW[tx_line, prd]), - value(m.Tx_Max_Capacity_MW[tx_line, prd]), + _finite_or_none(value(m.Tx_Min_Capacity_MW[tx_line, prd])), + _finite_or_none(value(m.Tx_Max_Capacity_MW[tx_line, prd])), ] for (tx_line, prd) in m.TX_OPR_PRDS ] diff --git a/gridpath/transmission/capacity/capacity_types/tx_spec.py b/gridpath/transmission/capacity/capacity_types/tx_spec.py index d6953129a0..d5632c1767 100644 --- a/gridpath/transmission/capacity/capacity_types/tx_spec.py +++ b/gridpath/transmission/capacity/capacity_types/tx_spec.py @@ -28,6 +28,7 @@ import csv import os.path +import pandas as pd from statistics import mean from pyomo.environ import Set, Param, Reals, NonNegativeReals @@ -42,10 +43,17 @@ write_validation_to_database, validate_dtypes, validate_idxs, - validate_missing_inputs, validate_column_monotonicity, ) +# A specified transmission line whose min (max) capacity is left blank in the +# inputs is treated as having no lower (upper) flow limit. The capacity params +# default to these sentinels, and the operational types skip the corresponding +# flow-limit constraint when the capacity is infinite (see +# min/max_limit_is_unconstrained_rule below). +Negative_Infinity = float("-inf") +Infinity = float("inf") + def add_model_components( m, @@ -116,8 +124,13 @@ def add_model_components( # Required Params ########################################################################### - m.tx_spec_min_cap_mw = Param(m.TX_SPEC_OPR_PRDS, within=Reals) - m.tx_spec_max_cap_mw = Param(m.TX_SPEC_OPR_PRDS, within=Reals) + # Optional caps: a blank min (max) in the inputs leaves the param at + # -inf (+inf), which the operational type reads as "no lower (upper) + # flow limit" and skips the corresponding constraint. + m.tx_spec_min_cap_mw = Param( + m.TX_SPEC_OPR_PRDS, within=Reals, default=Negative_Infinity + ) + m.tx_spec_max_cap_mw = Param(m.TX_SPEC_OPR_PRDS, within=Reals, default=Infinity) m.tx_spec_fixed_cost_per_mw_yr = Param( m.TX_SPEC_OPR_PRDS, within=NonNegativeReals, default=0 ) @@ -140,12 +153,30 @@ def max_transmission_capacity_rule(mod, tx, p): return mod.tx_spec_max_cap_mw[tx, p] +def min_limit_is_unconstrained_rule(mod, tx, p): + """Whether this line-period has no lower flow limit (blank min in inputs).""" + return mod.tx_spec_min_cap_mw[tx, p] == Negative_Infinity + + +def max_limit_is_unconstrained_rule(mod, tx, p): + """Whether this line-period has no upper flow limit (blank max in inputs).""" + return mod.tx_spec_max_cap_mw[tx, p] == Infinity + + def fixed_cost_rule(mod, g, p): """ The fixed cost of Tx lines of the *tx_spec* capacity type is a pre-specified number equal to the average capacity times the per-mw fixed cost for each of the project's operational periods. + + A line with no flow limit (infinite min or max capacity) has no + meaningful capacity to cost, so its fixed cost is zero. """ + if ( + mod.tx_spec_min_cap_mw[g, p] == Negative_Infinity + or mod.tx_spec_max_cap_mw[g, p] == Infinity + ): + return 0 return ( mean([abs(mod.tx_spec_min_cap_mw[g, p]), abs(mod.tx_spec_max_cap_mw[g, p])]) * mod.tx_spec_fixed_cost_per_mw_yr[g, p] @@ -167,32 +198,50 @@ def load_model_data( subproblem, stage, ): - data_portal.load( - filename=os.path.join( - scenario_directory, - weather_iteration, - hydro_iteration, - availability_iteration, - subproblem, - stage, - "inputs", - "specified_transmission_line_capacities.tab", - ), - select=( - "transmission_line", - "period", - "specified_tx_min_mw", - "specified_tx_max_mw", - "fixed_cost_per_mw_yr", - ), - index=m.TX_SPEC_OPR_PRDS, - param=( - m.tx_spec_min_cap_mw, - m.tx_spec_max_cap_mw, - m.tx_spec_fixed_cost_per_mw_yr, - ), + # min and max capacities are optional (a blank cell means "no flow limit + # in that direction"), so we cannot use a single data_portal.load() that + # ties index membership to parsing every param column. Instead we read the + # file manually, build TX_SPEC_OPR_PRDS from *every* row, and populate the + # capacity params per-cell, skipping blanks so they fall back to the + # ±Infinity defaults. This mirrors + # transmission/operations/transmission_flow_limits.py. + capacities_file = os.path.join( + scenario_directory, + weather_iteration, + hydro_iteration, + availability_iteration, + subproblem, + stage, + "inputs", + "specified_transmission_line_capacities.tab", ) + df = pd.read_csv(capacities_file, sep="\t") + + opr_prds = [] + min_cap = {} + max_cap = {} + fixed_cost = {} + for _, row in df.iterrows(): + tx = row["transmission_line"] + prd = int(row["period"]) + opr_prds.append((tx, prd)) + # "." (or a blank read as NaN) leaves the param at its ±inf default. + min_val = row["specified_tx_min_mw"] + if str(min_val) != "." and pd.notna(min_val): + min_cap[(tx, prd)] = float(min_val) + max_val = row["specified_tx_max_mw"] + if str(max_val) != "." and pd.notna(max_val): + max_cap[(tx, prd)] = float(max_val) + fc_val = row["fixed_cost_per_mw_yr"] + if str(fc_val) != "." and pd.notna(fc_val): + fixed_cost[(tx, prd)] = float(fc_val) + + data_portal.data()["TX_SPEC_OPR_PRDS"] = {None: opr_prds} + data_portal.data()["tx_spec_min_cap_mw"] = min_cap + data_portal.data()["tx_spec_max_cap_mw"] = max_cap + data_portal.data()["tx_spec_fixed_cost_per_mw_yr"] = fixed_cost + # Database ############################################################################### @@ -390,23 +439,14 @@ def validate_inputs( ), ) - # Check for missing values (vs. missing row entries above) - cols = ["min_mw", "max_mw"] - write_validation_to_database( - conn=conn, - scenario_id=scenario_id, - weather_iteration=weather_iteration, - hydro_iteration=hydro_iteration, - availability_iteration=availability_iteration, - subproblem_id=subproblem, - stage_id=stage, - gridpath_module=__name__, - db_table="inputs_transmission_specified_capacity", - severity="High", - errors=validate_missing_inputs(df, cols), - ) + # Note: min_mw and max_mw are intentionally NOT checked for missing + # values here -- a blank in either column is a valid input meaning "no + # flow limit in that direction" (the capacity param falls back to its + # ±Infinity default). - # check that min <= max + # check that min <= max (validate_column_monotonicity drops NaN rows, so + # lines with a blank/unconstrained min or max are skipped) + cols = ["min_mw", "max_mw"] write_validation_to_database( conn=conn, scenario_id=scenario_id, diff --git a/gridpath/transmission/operations/operational_types/tx_dcopf.py b/gridpath/transmission/operations/operational_types/tx_dcopf.py index 5f4ad7763e..ef08c47550 100644 --- a/gridpath/transmission/operations/operational_types/tx_dcopf.py +++ b/gridpath/transmission/operations/operational_types/tx_dcopf.py @@ -85,6 +85,21 @@ def add_model_components( | Two-dimensional set with transmission lines of the :code:`tx_dcopf` | | operational type and their operational timepoints. | +-------------------------------------------------------------------------+ + | | :code:`TX_DCOPF_OPR_TMPS_W_MIN_CONSTRAINT` | + | | + | Subset of :code:`TX_DCOPF_OPR_TMPS` restricted to line-timepoints whose | + | line-period has a lower flow limit (i.e. is in the transmission | + | capacity module's :code:`TX_OPR_PRDS_W_MIN_LIMIT`). The minimum-flow | + | constraint is built over this subset, so a line left unconstrained by | + | its capacity type gets no such constraint. | + +-------------------------------------------------------------------------+ + | | :code:`TX_DCOPF_OPR_TMPS_W_MAX_CONSTRAINT` | + | | + | Subset of :code:`TX_DCOPF_OPR_TMPS` restricted to line-timepoints whose | + | line-period has an upper flow limit (analogous to | + | :code:`TX_DCOPF_OPR_TMPS_W_MIN_CONSTRAINT`); scopes the maximum-flow | + | constraint. | + +-------------------------------------------------------------------------+ | @@ -211,6 +226,30 @@ def add_model_components( ), ) + # Operational timepoints whose line-period has a lower / upper flow limit; + # lines left unconstrained by their capacity type are excluded so no + # min/max flow constraint is built for them (see the transmission + # capacity module's TX_OPR_PRDS_W_MIN_LIMIT / _W_MAX_LIMIT). + m.TX_DCOPF_OPR_TMPS_W_MIN_CONSTRAINT = Set( + dimen=2, + within=m.TX_DCOPF_OPR_TMPS, + initialize=lambda mod: [ + (tx, tmp) + for (tx, tmp) in mod.TX_DCOPF_OPR_TMPS + if (tx, mod.period[tmp]) in mod.TX_OPR_PRDS_W_MIN_LIMIT + ], + ) + + m.TX_DCOPF_OPR_TMPS_W_MAX_CONSTRAINT = Set( + dimen=2, + within=m.TX_DCOPF_OPR_TMPS, + initialize=lambda mod: [ + (tx, tmp) + for (tx, tmp) in mod.TX_DCOPF_OPR_TMPS + if (tx, mod.period[tmp]) in mod.TX_OPR_PRDS_W_MAX_LIMIT + ], + ) + # Derived Sets ########################################################################### @@ -263,11 +302,11 @@ def add_model_components( ########################################################################### m.TxDcopf_Min_Transmit_Constraint = Constraint( - m.TX_DCOPF_OPR_TMPS, rule=min_transmit_rule + m.TX_DCOPF_OPR_TMPS_W_MIN_CONSTRAINT, rule=min_transmit_rule ) m.TxDcopf_Max_Transmit_Constraint = Constraint( - m.TX_DCOPF_OPR_TMPS, rule=max_transmit_rule + m.TX_DCOPF_OPR_TMPS_W_MAX_CONSTRAINT, rule=max_transmit_rule ) m.TxDcopf_Kirchhoff_Voltage_Law_Constraint = Constraint( diff --git a/gridpath/transmission/operations/operational_types/tx_simple.py b/gridpath/transmission/operations/operational_types/tx_simple.py index 88bd7e8a2f..14d2e671c8 100644 --- a/gridpath/transmission/operations/operational_types/tx_simple.py +++ b/gridpath/transmission/operations/operational_types/tx_simple.py @@ -68,6 +68,21 @@ def add_model_components( | Two-dimensional set with transmission lines of the :code:`tx_simple` | | operational type and their operational timepoints. | +-------------------------------------------------------------------------+ + | | :code:`TX_SIMPLE_OPR_TMPS_W_MIN_LIMIT` | + | | + | Subset of :code:`TX_SIMPLE_OPR_TMPS` restricted to line-timepoints | + | whose line-period has a lower flow limit (i.e. is in the transmission | + | capacity module's :code:`TX_OPR_PRDS_W_MIN_LIMIT`). The minimum-flow | + | and "from"-direction loss constraints are built over this subset, so a | + | line left unconstrained by its capacity type gets no such constraint. | + +-------------------------------------------------------------------------+ + | | :code:`TX_SIMPLE_OPR_TMPS_W_MAX_LIMIT` | + | | + | Subset of :code:`TX_SIMPLE_OPR_TMPS` restricted to line-timepoints | + | whose line-period has an upper flow limit (analogous to | + | :code:`TX_SIMPLE_OPR_TMPS_W_MIN_LIMIT`); scopes the maximum-flow and | + | "to"-direction loss constraints. | + +-------------------------------------------------------------------------+ +-------------------------------------------------------------------------+ | Params | @@ -180,6 +195,36 @@ def add_model_components( ), ) + # Operational timepoints whose line-period has a lower / upper flow limit. + # Lines left unconstrained by their capacity type (e.g. a tx_spec line with + # a blank min or max) are excluded, so no min/max flow constraint is built + # for them. TX_OPR_PRDS_W_MIN_LIMIT / _W_MAX_LIMIT come from the + # transmission capacity module. + # Note: distinct from the identically-purposed but separately-fed + # TX_SIMPLE_OPR_TMPS_W_{MIN,MAX}_CONSTRAINT sets in + # transmission/operations/transmission_flow_limits.py (which come from the + # optional transmission_flow_limits inputs). These "_LIMIT" sets come from + # the line's *capacity* and gate the capacity-based transmit constraints. + m.TX_SIMPLE_OPR_TMPS_W_MIN_LIMIT = Set( + dimen=2, + within=m.TX_SIMPLE_OPR_TMPS, + initialize=lambda mod: [ + (tx, tmp) + for (tx, tmp) in mod.TX_SIMPLE_OPR_TMPS + if (tx, mod.period[tmp]) in mod.TX_OPR_PRDS_W_MIN_LIMIT + ], + ) + + m.TX_SIMPLE_OPR_TMPS_W_MAX_LIMIT = Set( + dimen=2, + within=m.TX_SIMPLE_OPR_TMPS, + initialize=lambda mod: [ + (tx, tmp) + for (tx, tmp) in mod.TX_SIMPLE_OPR_TMPS + if (tx, mod.period[tmp]) in mod.TX_OPR_PRDS_W_MAX_LIMIT + ], + ) + # Params ########################################################################### m.tx_simple_loss_factor = Param(m.TX_SIMPLE, within=PercentFraction, default=0) @@ -196,11 +241,11 @@ def add_model_components( ########################################################################### m.TxSimple_Min_Transmit_Constraint = Constraint( - m.TX_SIMPLE_OPR_TMPS, rule=min_transmit_rule + m.TX_SIMPLE_OPR_TMPS_W_MIN_LIMIT, rule=min_transmit_rule ) m.TxSimple_Max_Transmit_Constraint = Constraint( - m.TX_SIMPLE_OPR_TMPS, rule=max_transmit_rule + m.TX_SIMPLE_OPR_TMPS_W_MAX_LIMIT, rule=max_transmit_rule ) m.TxSimple_Losses_LZ_From_Constraint = Constraint( @@ -211,12 +256,15 @@ def add_model_components( m.TX_SIMPLE_OPR_TMPS, rule=losses_lz_to_rule ) + # The loss upper bounds are the flow capacity times the loss factor, so + # they only apply where that capacity is finite (min for the "from" + # direction, max for the "to" direction). m.TxSimple_Max_Losses_From_Constraint = Constraint( - m.TX_SIMPLE_OPR_TMPS, rule=max_losses_from_rule + m.TX_SIMPLE_OPR_TMPS_W_MIN_LIMIT, rule=max_losses_from_rule ) m.TxSimple_Max_Losses_To_Constraint = Constraint( - m.TX_SIMPLE_OPR_TMPS, rule=max_losses_to_rule + m.TX_SIMPLE_OPR_TMPS_W_MAX_LIMIT, rule=max_losses_to_rule ) diff --git a/gridpath/transmission/operations/operational_types/tx_simple_binary.py b/gridpath/transmission/operations/operational_types/tx_simple_binary.py index 3163c6163e..31a328d749 100644 --- a/gridpath/transmission/operations/operational_types/tx_simple_binary.py +++ b/gridpath/transmission/operations/operational_types/tx_simple_binary.py @@ -70,19 +70,22 @@ def add_model_components( | Two-dimensional set with transmission lines of the :code:`tx_simple_binary` | | operational type and their operational timepoints. | +-------------------------------------------------------------------------+ - | | :code:`TX_SIMPLE_BINARY_OPR_TMPS_W_MIN_CONSTRAINT` | + | | :code:`TX_SIMPLE_BINARY_OPR_TMPS_W_MIN_CONSTRAINT` | | | - | Two-dimensional set with transmission lines of the :code:`tx_simple_binary` | - | operational type and their operational timepoints to describe all | - | possible transmission-timepoint combinations for transmission lines | - | with a minimum flow specified. | + | Subset of :code:`TX_SIMPLE_BINARY_OPR_TMPS` restricted to | + | line-timepoints whose line-period has a lower flow limit (i.e. is in | + | the transmission capacity module's :code:`TX_OPR_PRDS_W_MIN_LIMIT`). | + | The minimum-flow, negative-direction big-M, and "from"-direction loss | + | constraints are built over this subset, so a line left unconstrained by | + | its capacity type gets no such constraint. | +-------------------------------------------------------------------------+ - | | :code:`TX_SIMPLE_BINARY_OPR_TMPS_W_MAX_CONSTRAINT` | + | | :code:`TX_SIMPLE_BINARY_OPR_TMPS_W_MAX_CONSTRAINT` | | | - | Two-dimensional set with transmission lines of the :code:`tx_simple_binary` | - | operational type and their operational timepoints to describe all | - | possible transmission-timepoint combinations for transmission lines | - | with a maximum flow specified. | + | Subset of :code:`TX_SIMPLE_BINARY_OPR_TMPS` restricted to | + | line-timepoints whose line-period has an upper flow limit (analogous to | + | :code:`TX_SIMPLE_BINARY_OPR_TMPS_W_MIN_CONSTRAINT`); scopes the | + | maximum-flow, positive-direction big-M, and "to"-direction loss | + | constraints. | +-------------------------------------------------------------------------+ +-------------------------------------------------------------------------+ @@ -199,12 +202,28 @@ def add_model_components( ), ) + # Operational timepoints whose line-period has a lower / upper flow limit; + # lines left unconstrained by their capacity type are excluded so no + # min/max (or directional big-M) constraint is built for them (see the + # transmission capacity module's TX_OPR_PRDS_W_MIN_LIMIT / _W_MAX_LIMIT). m.TX_SIMPLE_BINARY_OPR_TMPS_W_MIN_CONSTRAINT = Set( - dimen=2, within=m.TX_SIMPLE_BINARY_OPR_TMPS + dimen=2, + within=m.TX_SIMPLE_BINARY_OPR_TMPS, + initialize=lambda mod: [ + (tx, tmp) + for (tx, tmp) in mod.TX_SIMPLE_BINARY_OPR_TMPS + if (tx, mod.period[tmp]) in mod.TX_OPR_PRDS_W_MIN_LIMIT + ], ) m.TX_SIMPLE_BINARY_OPR_TMPS_W_MAX_CONSTRAINT = Set( - dimen=2, within=m.TX_SIMPLE_BINARY_OPR_TMPS + dimen=2, + within=m.TX_SIMPLE_BINARY_OPR_TMPS, + initialize=lambda mod: [ + (tx, tmp) + for (tx, tmp) in mod.TX_SIMPLE_BINARY_OPR_TMPS + if (tx, mod.period[tmp]) in mod.TX_OPR_PRDS_W_MAX_LIMIT + ], ) # Params @@ -248,20 +267,27 @@ def binary_transmit_power_rule(mod, tx, tmp): # Constraints ########################################################################### + # The directional big-M constraints use the flow capacity as the big-M + # (binary * capacity), so they only apply where that capacity is finite: + # positive direction uses the max capacity, negative uses the min. A line + # left unconstrained in a direction has no big-M to enforce, so its + # directional constraint is skipped (the binary cannot prevent simultaneous + # bidirectional flow on a limitless line -- an inherent, documented + # limitation of pairing tx_simple_binary with an unconstrained line). m.TxSimpleBinary_Positive_Direction_Constraint = Constraint( - m.TX_SIMPLE_BINARY_OPR_TMPS, rule=positive_direction_rule + m.TX_SIMPLE_BINARY_OPR_TMPS_W_MAX_CONSTRAINT, rule=positive_direction_rule ) m.TxSimpleBinary_Negative_Direction_Constraint = Constraint( - m.TX_SIMPLE_BINARY_OPR_TMPS, rule=negative_direction_rule + m.TX_SIMPLE_BINARY_OPR_TMPS_W_MIN_CONSTRAINT, rule=negative_direction_rule ) m.TxSimpleBinary_Min_Transmit_Constraint = Constraint( - m.TX_SIMPLE_BINARY_OPR_TMPS, rule=min_transmit_rule + m.TX_SIMPLE_BINARY_OPR_TMPS_W_MIN_CONSTRAINT, rule=min_transmit_rule ) m.TxSimpleBinary_Max_Transmit_Constraint = Constraint( - m.TX_SIMPLE_BINARY_OPR_TMPS, rule=max_transmit_rule + m.TX_SIMPLE_BINARY_OPR_TMPS_W_MAX_CONSTRAINT, rule=max_transmit_rule ) m.TxSimpleBinary_Losses_LZ_From_Constraint = Constraint( @@ -272,12 +298,14 @@ def binary_transmit_power_rule(mod, tx, tmp): m.TX_SIMPLE_BINARY_OPR_TMPS, rule=losses_lz_to_rule ) + # Loss upper bounds are the flow capacity times the loss factor, so they + # only apply where that capacity is finite (min for "from", max for "to"). m.TxSimpleBinary_Max_Losses_From_Constraint = Constraint( - m.TX_SIMPLE_BINARY_OPR_TMPS, rule=max_losses_from_rule + m.TX_SIMPLE_BINARY_OPR_TMPS_W_MIN_CONSTRAINT, rule=max_losses_from_rule ) m.TxSimpleBinary_Max_Losses_To_Constraint = Constraint( - m.TX_SIMPLE_BINARY_OPR_TMPS, rule=max_losses_to_rule + m.TX_SIMPLE_BINARY_OPR_TMPS_W_MAX_CONSTRAINT, rule=max_losses_to_rule ) diff --git a/tests/test_scaling.py b/tests/test_scaling.py new file mode 100644 index 0000000000..3091f5fe61 --- /dev/null +++ b/tests/test_scaling.py @@ -0,0 +1,297 @@ +# 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. + +""" +Tests for numerical scaling (gridpath/auxiliary/scaling.py) and its integration +into the solve path (gridpath/run_scenario.py::solve_problem). + +The correctness property being verified: solving with scale factors and mapping +the solution back yields the same objective, variable values, and duals as +solving unscaled -- i.e. the scaling is an exact reformulation, not an +approximation. +""" + +import tempfile +import types +import unittest + +from pyomo.environ import ( + Binary, + ConcreteModel, + Constraint, + NonNegativeReals, + Objective, + Reals, + Set, + SolverFactory, + Suffix, + TransformationFactory, + Var, + minimize, + value, +) + +from gridpath import run_scenario +from gridpath.auxiliary.scaling import ( + assign_scaling_factors, + classify_variable_units, +) + +# Scale factors used throughout: 1000 (MW -> GW) and 1e6 ($ -> $M). +POWER_SCALE = 1000.0 +DOLLAR_SCALE = 1e6 + + +def _build_dispatch_model(): + """Build a small LP that mimics GridPath's structure and naming. + + A two-timepoint economic dispatch: a generator (power, MW) and unserved + energy (power, MW) serve a load; a cost variable (dollars) is defined by a + $/MWh cost curve. The load-balance constraint is binding, so its dual is a + non-zero marginal price -- exercising dual propagation through the + S_dollar/S_power ratio. + + Returns: + A Pyomo ``ConcreteModel`` with a ``dual`` Suffix (as run_scenario sets). + """ + m = ConcreteModel() + m.T = Set(initialize=[1, 2], ordered=True) + load = {1: 800.0, 2: 5000.0} # tmp 2 exceeds the 4000 MW gen cap -> unserved + + m.GenSimple_Provide_Power_MW = Var(m.T, bounds=(0, 4000)) + m.Unserved_Energy_MW = Var(m.T, within=NonNegativeReals) + m.Variable_OM_Curve_Cost = Var(m.T, within=Reals) + + def meet_load_rule(mod, t): + return mod.GenSimple_Provide_Power_MW[t] + mod.Unserved_Energy_MW[t] == load[t] + + m.Meet_Load_Constraint = Constraint(m.T, rule=meet_load_rule) + + # Cost curve: $30/MWh generation, $500/MWh unserved-energy penalty. The + # rate constants live inside this (cost-definition) row. + def cost_rule(mod, t): + return mod.Variable_OM_Curve_Cost[t] >= ( + 30 * mod.GenSimple_Provide_Power_MW[t] + 500 * mod.Unserved_Energy_MW[t] + ) + + m.Cost_Constraint = Constraint(m.T, rule=cost_rule) + + m.NPV = Objective( + expr=sum(m.Variable_OM_Curve_Cost[t] for t in m.T), sense=minimize + ) + # run_scenario declares this on every model; scale_model uses it to map + # duals back. + m.dual = Suffix(direction=Suffix.IMPORT) + return m + + +def _fake_parsed_arguments( + power_scale_factor, dollar_scale_factor, scale_mode="out_of_place" +): + """Build a minimal parsed-arguments stand-in for solve_problem/solve. + + Points the (non-existent) scenario directory at a temp location so that no + ``solver_options.csv`` is found and ``solve`` falls back to cbc. + """ + return types.SimpleNamespace( + quiet=True, + power_scale_factor=power_scale_factor, + dollar_scale_factor=dollar_scale_factor, + scale_mode=scale_mode, + solver=None, + solver_executable=None, + mute_solver_output=True, + keepfiles=False, + symbolic=False, + scenario_location=tempfile.gettempdir(), + scenario="_scaling_unittest_no_such_scenario", + ) + + +class TestClassifyVariableUnits(unittest.TestCase): + """Name-based unit classification against real GridPath variable names.""" + + def test_power_energy_names(self): + for name in [ + "GenSimple_Provide_Power_MW", + "Stor_Starting_Energy_in_Storage_MWh", + "Inertia_Reserves_Violation_MWs", + "Transmission_Target_Energy_MW_Neg_Dir", + "Net_Market_Purchased_Power", # unsuffixed, trailing "Power" + ]: + self.assertEqual(classify_variable_units(name), "power", msg=name) + + def test_dollar_names(self): + for name in [ + "Hurdle_Cost_Pos_Dir", # "Cost" token, not last + "Variable_OM_Curve_Cost", + "Carbon_Tax_Cost", + "Ramp_Up_Tuning_Cost", + ]: + self.assertEqual(classify_variable_units(name), "dollar", msg=name) + + def test_unrecognized_names(self): + for name in [ + "Fuel_Prod_Consume_Power_PowerUnit", # trailing "PowerUnit", not "Power" + "Import_Carbon_Emissions_Tons", + "Load_Component_Modifier_Fraction_Invested", + "LZ_Exports", + ]: + self.assertIsNone(classify_variable_units(name), msg=name) + + +class TestAssignScalingFactors(unittest.TestCase): + """Factor assignment on a built instance.""" + + def setUp(self): + self.m = _build_dispatch_model() + assign_scaling_factors(self.m, POWER_SCALE, DOLLAR_SCALE) + + def test_positive_factor_validation(self): + for bad in [(0, 1), (1, 0), (-1, 1), (1, -5)]: + with self.assertRaises(ValueError): + assign_scaling_factors(_build_dispatch_model(), bad[0], bad[1]) + + def test_power_variable_factor(self): + self.assertEqual( + self.m.scaling_factor[self.m.GenSimple_Provide_Power_MW], + 1.0 / POWER_SCALE, + ) + self.assertEqual( + self.m.scaling_factor[self.m.Unserved_Energy_MW], 1.0 / POWER_SCALE + ) + + def test_dollar_variable_factor(self): + self.assertEqual( + self.m.scaling_factor[self.m.Variable_OM_Curve_Cost], + 1.0 / DOLLAR_SCALE, + ) + + def test_objective_factor(self): + self.assertEqual(self.m.scaling_factor[self.m.NPV], 1.0 / DOLLAR_SCALE) + + def test_homogeneous_power_constraint_factor(self): + # Meet_Load: sum(MW) == load -> power factor. + self.assertEqual( + self.m.scaling_factor[self.m.Meet_Load_Constraint], + 1.0 / POWER_SCALE, + ) + + def test_cost_definition_constraint_factor(self): + # Cost_Constraint contains a dollar variable -> dollar factor. + self.assertEqual( + self.m.scaling_factor[self.m.Cost_Constraint], 1.0 / DOLLAR_SCALE + ) + + def test_container_granularity(self): + # One suffix entry per container (not per index), so the count is small + # and independent of the timepoint set size: 3 vars + 2 constraints + 1 + # objective = 6 entries, even though each indexed component has 2 data. + self.assertEqual(len(self.m.scaling_factor), 6) + + def test_integer_variable_not_scaled(self): + m = _build_dispatch_model() + m.GenNewBin_Build = Var(within=Binary) + assign_scaling_factors(m, POWER_SCALE, DOLLAR_SCALE) + self.assertNotIn(m.GenNewBin_Build, m.scaling_factor) + + +class TestScaledSolveEquivalence(unittest.TestCase): + """The scaled solve must reproduce the unscaled solution exactly.""" + + @classmethod + def setUpClass(cls): + if not SolverFactory("cbc").available(): + raise unittest.SkipTest("cbc not available") + + def _solve( + self, power_scale_factor, dollar_scale_factor, scale_mode="out_of_place" + ): + instance = _build_dispatch_model() + args = _fake_parsed_arguments( + power_scale_factor, dollar_scale_factor, scale_mode=scale_mode + ) + solved, _ = run_scenario.solve_problem(args, instance) + gen = {t: value(solved.GenSimple_Provide_Power_MW[t]) for t in solved.T} + use = {t: value(solved.Unserved_Energy_MW[t]) for t in solved.T} + price = {t: solved.dual[solved.Meet_Load_Constraint[t]] for t in solved.T} + return value(solved.NPV), gen, use, price + + def _check_matches_native(self, base, scaled): + base_obj, base_gen, base_use, base_price = base + scl_obj, scl_gen, scl_use, scl_price = scaled + self.assertAlmostEqual(base_obj, scl_obj, places=3) + for t in base_gen: + self.assertAlmostEqual(base_gen[t], scl_gen[t], places=4, msg=f"gen[{t}]") + self.assertAlmostEqual(base_use[t], scl_use[t], places=4, msg=f"use[{t}]") + # Marginal price ($/MWh) recovered via the S_dollar/S_power ratio. + self.assertAlmostEqual( + base_price[t], scl_price[t], places=4, msg=f"price[{t}]" + ) + + def test_out_of_place_matches_native(self): + base = self._solve(1.0, 1.0) + scaled = self._solve(POWER_SCALE, DOLLAR_SCALE, scale_mode="out_of_place") + # Objective (native dollars): 800*30 + 4000*30 + 1000*500 = 644000. + self.assertAlmostEqual(base[0], 644000.0, places=1) + self._check_matches_native(base, scaled) + # tmp 2 is short 1000 MW -> its marginal price is the $500/MWh penalty. + self.assertAlmostEqual(scaled[3][2], 500.0, places=3) + + def test_in_place_matches_native(self): + base = self._solve(1.0, 1.0) + scaled = self._solve(POWER_SCALE, DOLLAR_SCALE, scale_mode="in_place") + self._check_matches_native(base, scaled) + self.assertAlmostEqual(scaled[3][2], 500.0, places=3) + + def test_no_scaling_path_leaves_no_suffix(self): + # The (1.0, 1.0) default path must not touch the instance (no clone, no + # scaling_factor suffix). + instance = _build_dispatch_model() + args = _fake_parsed_arguments(1.0, 1.0) + solved, _ = run_scenario.solve_problem(args, instance) + self.assertIs(solved, instance) + self.assertFalse(hasattr(instance, "scaling_factor")) + + +class TestIncompatibleFlagGuard(unittest.TestCase): + """Scaling combined with a load-solution / lp-only flag must fail fast.""" + + def test_scaling_with_load_solution_raises(self): + # The guard runs right after argument parsing, before the scenario + # directory is checked, so a placeholder --scenario is enough. + with self.assertRaises(ValueError): + run_scenario.main( + [ + "--scenario", + "_scaling_unittest", + "--power_scale_factor", + "1000", + "--load_highs_solution", + ] + ) + + def test_scaling_with_lp_only_raises(self): + with self.assertRaises(ValueError): + run_scenario.main( + [ + "--scenario", + "_scaling_unittest", + "--dollar_scale_factor", + "1000000", + "--create_lp_problem_file_only", + ] + ) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/test_scaling_examples.py b/tests/test_scaling_examples.py new file mode 100644 index 0000000000..804f1dfd6a --- /dev/null +++ b/tests/test_scaling_examples.py @@ -0,0 +1,188 @@ +# 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. + +""" +End-to-end equivalence test for numerical scaling, built on the example suite. + +``TestScaledExamples`` subclasses ``tests.test_examples.TestExamples`` and +re-runs every example scenario with the power/energy and dollar scale factors +turned on, checking that the objective still matches the same stored expected +value (to a relative tolerance). Because it only overrides the comparison hook +(``run_and_check_objective``), it automatically inherits the parent's test +database setup and *every* ``test_*`` method -- including any example test added +later. A new scenario that happens to break under scaling therefore fails here +without anyone having to remember to add it. + +Rationale for the relative tolerance: the unscaled suite asserts an *absolute* +objective match (``places=1``), which only holds because GridPath is +bit-deterministic per platform. Scaling deliberately changes the floating-point +numbers the solver sees, so identical answers can differ in the last few digits +(and, for degenerate LPs, the solver may settle on a different equally-optimal +vertex, differing by ~1e-7 relative). A relative tolerance passes those while +still catching real breakage, which shows up as a large relative difference (or +an outright error, as an unhandled missing-dual once did for ``test_markets``). + +This complements the fast, DB-free unit tests in ``test_scaling.py`` (classifier, +factor assignment, exact round-trip, argument guards), which pinpoint *why* +something breaks; this suite is the broad integration net. +""" + +import multiprocessing + +# Import the module (not the TestExamples class by name): pytest collects any +# TestCase subclass bound at module scope, so importing the class directly would +# make the unscaled TestExamples suite run a second time in this file. Referencing +# it through the module keeps only TestScaledExamples collected here. +from tests import test_examples +from tests.test_examples import ( + objective_function_overwrite, + DB_PATH, + EXAMPLES_DIRECTORY, +) +from gridpath import run_end_to_end + + +class TestScaledExamples(test_examples.TestExamples): + """Run the full example suite with numerical scaling on.""" + + # MW -> GW, MWh -> GWh; $ -> $M. + POWER_SCALE_FACTOR = 1000.0 + DOLLAR_SCALE_FACTOR = 1000000.0 + + # Which scaling implementation to exercise ("out_of_place" or "in_place"). + # Both must give identical native-unit results; the unit tests in + # test_scaling.py check that equivalence directly. This suite runs the + # default (out_of_place) mode. To sweep the full example set through the + # in_place path, subclass and flip this attribute: + # class TestScaledExamplesInPlace(TestScaledExamples): + # SCALE_MODE = "in_place" + # (not added as an always-on suite -- it would roughly triple example + # runtime for little marginal signal over the unit-level equivalence tests). + SCALE_MODE = "out_of_place" + + # Relative tolerance on the objective: |actual - expected| <= tol * |expected|. + # Most scenarios match to machine precision (~1e-13). Degenerate LPs, though, + # can settle on a different equally-optimal vertex once scaling perturbs the + # numbers -- especially scenarios dominated by commodity chains the classifier + # leaves unscaled (fuel, emissions, water), which get only partial + # conditioning. The observed worst case across the example suite is ~2.85e-6 + # (test_new_solar_carbon_cap_2zones_tx_hydrogen_prod_new); 1e-5 covers that + # with margin while real breakage still fails loudly (a crash, or a gross + # objective difference). + RELATIVE_TOLERANCE = 1e-5 + # Absolute floor for the (not expected in practice) zero-objective case. + ABSOLUTE_FLOOR = 1e-6 + + def check_validation(self, test): + """Skip input validation in the scaled suite. + + Validation runs against the same inputs regardless of solve-time + scaling, so it is already covered by ``TestExamples``; re-running it here + would only double the work. + """ + pass + + def _assert_objective_close(self, expected, actual, msg=None): + """Recursively compare (possibly nested per-subproblem/stage) objectives + by relative tolerance.""" + if isinstance(expected, dict): + self.assertIsInstance(actual, dict, msg) + self.assertEqual(expected.keys(), actual.keys(), msg) + for key in expected: + self._assert_objective_close(expected[key], actual[key], msg=msg) + else: + bound = ( + self.RELATIVE_TOLERANCE * max(abs(expected), abs(actual)) + + self.ABSOLUTE_FLOOR + ) + self.assertLessEqual( + abs(expected - actual), + bound, + msg=( + f"{msg or ''} expected={expected} actual={actual} " + f"rel_diff={abs(expected - actual) / (abs(expected) or 1):.3e} " + f"(tol={self.RELATIVE_TOLERANCE})" + ), + ) + + def run_and_check_objective( + self, + scenario_name, + expected_objective, + additional_args=[], + solver=None, + parallel=1, + ): + """Run a scenario with scaling on and check the objective by relative + tolerance. + + Mirrors ``TestExamples.run_and_check_objective`` but (a) appends the + scale-factor arguments, (b) compares with a relative tolerance, and (c) + does not write the ``actual_objective`` column back to the tracked + expected-values CSV (that column is for the unscaled baseline). + """ + args_to_pass = [ + "--database", + DB_PATH, + "--scenario", + scenario_name, + "--scenario_location", + EXAMPLES_DIRECTORY, + "--n_parallel_get_inputs", + str(parallel), + "--n_parallel_solve", + str(parallel), + "--quiet", + "--mute_solver_output", + "--testing", + "--power_scale_factor", + str(self.POWER_SCALE_FACTOR), + "--dollar_scale_factor", + str(self.DOLLAR_SCALE_FACTOR), + "--scale_mode", + self.SCALE_MODE, + ] + additional_args + + if solver is not None: + args_to_pass.append("--solver") + args_to_pass.append(solver) + + actual_objective = run_end_to_end.main(args_to_pass) + + expected_objective = objective_function_overwrite( + scenario_name=scenario_name, starting_objective=expected_objective + ) + + # Flatten a multiprocessing manager proxy dict to a plain dict (only + # relevant when parallel > 1); mirrors the parent. + if hasattr(multiprocessing, "managers"): + if isinstance(actual_objective, multiprocessing.managers.DictProxy): + actual_objective_copy = dict(actual_objective.copy()) + for subproblem in actual_objective.keys(): + if isinstance( + actual_objective[subproblem], + multiprocessing.managers.DictProxy, + ): + actual_objective_copy[subproblem] = dict( + actual_objective_copy[subproblem].copy() + ) + actual_objective = actual_objective_copy + + self._assert_objective_close( + expected_objective, actual_objective, msg=scenario_name + ) + + +if __name__ == "__main__": + import unittest + + unittest.main() diff --git a/tests/transmission/capacity/capacity_types/test_tx_spec.py b/tests/transmission/capacity/capacity_types/test_tx_spec.py index 9b41f6ffe5..c224f0d335 100644 --- a/tests/transmission/capacity/capacity_types/test_tx_spec.py +++ b/tests/transmission/capacity/capacity_types/test_tx_spec.py @@ -16,7 +16,9 @@ from collections import OrderedDict from importlib import import_module import os.path +import shutil import sys +import tempfile import unittest from tests.common_functions import create_abstract_model, add_components_and_load_data @@ -191,6 +193,59 @@ def test_data_loaded_correctly(self): ) self.assertDictEqual(expected_fixed_cost, actual_fixed_cost) + def test_blank_caps_are_unconstrained(self): + """A blank min/max leaves the line-period in TX_SPEC_OPR_PRDS but sets + the capacity param to its +/-inf default (i.e. no flow limit).""" + tmp_dir = tempfile.mkdtemp() + try: + staged = os.path.join(tmp_dir, "test_data") + shutil.copytree(TEST_DATA_DIRECTORY, staged) + cap_file = os.path.join( + staged, "inputs", "specified_transmission_line_capacities.tab" + ) + # Blank Tx1's max (2020) and Tx2's min (2020); leave the rows in place. + rows = open(cap_file).read().split("\n") + out = [] + for r in rows: + f = r.split("\t") + if f[0] == "Tx1" and f[1] == "2020": + f[3] = "." # specified_tx_max_mw -> no upper limit + r = "\t".join(f) + elif f[0] == "Tx2" and f[1] == "2020": + f[2] = "." # specified_tx_min_mw -> no lower limit + r = "\t".join(f) + out.append(r) + open(cap_file, "w").write("\n".join(out)) + + m, data = add_components_and_load_data( + prereq_modules=IMPORTED_PREREQ_MODULES, + module_to_test=MODULE_BEING_TESTED, + test_data_dir=staged, + weather_iteration="", + hydro_iteration="", + availability_iteration="", + subproblem="", + stage="", + ) + instance = m.create_instance(data) + + # Row still present (so the line stays operational / in TX_OPR_PRDS). + self.assertIn(("Tx1", 2020), list(instance.TX_SPEC_OPR_PRDS)) + self.assertIn(("Tx2", 2020), list(instance.TX_SPEC_OPR_PRDS)) + + # Blank cells fall back to the +/-inf defaults. + self.assertEqual( + instance.tx_spec_max_cap_mw["Tx1", 2020], float("inf") + ) + self.assertEqual( + instance.tx_spec_min_cap_mw["Tx2", 2020], float("-inf") + ) + # The non-blank direction on each line is unaffected. + self.assertEqual(instance.tx_spec_min_cap_mw["Tx1", 2020], -10) + self.assertEqual(instance.tx_spec_max_cap_mw["Tx2", 2020], 10) + finally: + shutil.rmtree(tmp_dir, ignore_errors=True) + if __name__ == "__main__": unittest.main()