Problem Scaling

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
 
 

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,

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]

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
     ]