Support for loading HiGHS solutions + tweaks to CPLEX solution parsing

Author: anamileva

gridpath/common_functions.py

@@ -1,4 +1,5 @@
-# Copyright 2016-2023 Blue Marble Analytics LLC.
+# Copyright 2016-2025 Blue Marble Analytics LLC.
+# Copyright 2026 Sylvan Energy Analytics LLC.
 #
 # Licensed under the Apache License, Version 2.0 (the "License");
 # you may not use this file except in compliance with the License.
@@ -205,6 +206,12 @@ def get_run_scenario_parser():
         action="store_true",
         help="Skip solve and load results from a Gurobi solution file instead.",
     )
+    parser.add_argument(
+        "--load_highs_solution",
+        default=False,
+        action="store_true",
+        help="Skip solve and load results from a HiGHS solution file instead.",
+    )
     # Solver options
     parser.add_argument(
         "--solver",

gridpath/run_scenario.py

@@ -1,4 +1,5 @@
-# Copyright 2016-2023 Blue Marble Analytics LLC.
+# Copyright 2016-2025 Blue Marble Analytics LLC.
+# Copyright 2026 Sylvan Energy Analytics LLC.
 #
 # Licensed under the Apache License, Version 2.0 (the "License");
 # you may not use this file except in compliance with the License.
@@ -320,6 +321,11 @@ def run_optimization_for_subproblem_stage(
                 prob_sol_files_directory=prob_sol_files_directory,
                 solution_filename="gurobi_solution.json",
             )
+        if parsed_arguments.load_highs_solution:
+            solved_instance, results, dynamic_components = load_highs_xml_solution(
+                prob_sol_files_directory=prob_sol_files_directory,
+                solution_filename="highs_solution.sol",
+            )
         else:
             dynamic_components, instance = create_problem(
                 scenario_directory=scenario_directory,
@@ -357,7 +363,7 @@ def run_optimization_for_subproblem_stage(
                 symbol_map = instance.solutions.symbol_map[smap_id]
 
                 symbol_cuid_pairs = tuple(
-                    (symbol, ComponentUID(var_weakref(), cuid_buffer={}))
+                    (symbol, ComponentUID(var_weakref, cuid_buffer={}))
                     for symbol, var_weakref in symbol_map.bySymbol.items()
                 )
 
@@ -1638,8 +1644,9 @@ def load_cplex_xml_solution(
             type_tag.get("index"),
             type_tag.get("value"),
         )
-        if not var_id == "ONE_VAR_CONSTANT":
-            symbol_map.bySymbol[var_id]().value = float(value)
+        # "x2" with value None added for CPLEXSolution version 1.2
+        if not var_id in ["ONE_VAR_CONSTANT", "x2"]:
+            symbol_map.bySymbol[var_id].value = float(value)
 
     # Constraints
     for type_tag in root.findall("linearConstraints/constraint"):
@@ -1649,8 +1656,9 @@ def load_cplex_xml_solution(
             type_tag.get("dual"),
         )
         if not constraint_id_w_extra_symbols == "c_e_ONE_VAR_CONSTANT":
-            constraint_id = constraint_id_w_extra_symbols[4:-1]
-            instance.dual[symbol_map.bySymbol[constraint_id]()] = float(dual)
+            # constraint_id = constraint_id_w_extra_symbols[4:-1]
+            constraint_id = constraint_id_w_extra_symbols
+            instance.dual[symbol_map.bySymbol[constraint_id]] = float(dual)
 
     # Solver status
     header = root.findall("header")[0]  # Need a check that there is only one element
@@ -1715,6 +1723,95 @@ def load_gurobi_json_solution(
     return instance, results, dynamic_components
 
 
+def load_highs_xml_solution(
+    prob_sol_files_directory, solution_filename="highs_solution.sol"
+):
+    """
+    :param prob_sol_files_directory:
+    :param solution_filename:
+    :return:
+    """
+    print(
+        "Loading results from solution file {}...".format(
+            os.path.join(prob_sol_files_directory, solution_filename)
+        )
+    )
+    instance, dynamic_components, symbol_map = load_problem_info(
+        prob_sol_files_directory=prob_sol_files_directory
+    )
+
+    # Read HiGHS solution file
+    with open(os.path.join(prob_sol_files_directory, solution_filename), "r") as f:
+        lines = f.readlines()
+
+    # Model status is  th second line
+    model_status = lines[1].strip() if len(lines) > 1 else "Unknown"
+    termination_condition = model_status.lower()
+    solver_status = "ok" if termination_condition == "optimal" else "unknown"
+
+    # Parse the HiGHS solution file
+    section = None
+
+    for line in lines:
+        line = line.strip()
+
+        # Skip empty lines and comments
+        if not line or line.startswith("#"):
+            # Check for section headers in comments
+            if line == "# Primal solution values":
+                section = "primal_start"
+            elif "# Columns" in line and section == "primal_start":
+                section = "primal_columns"
+            elif "# Rows" in line and section == "primal_columns":
+                section = "primal_rows"
+            elif line == "# Dual solution values":
+                section = "dual_start"
+            elif "# Columns" in line and section == "dual_start":
+                section = "dual_columns"
+            elif "# Rows" in line and section in ["dual_start", "dual_columns"]:
+                section = "dual_rows"
+            elif line == "# Basis":
+                # Stop parsing once we reach basis section
+                break
+            continue
+
+        # Skip non-data lines
+        if line in ["Model status", "Feasible", "Valid"]:
+            continue
+        if line.startswith("Objective "):
+            continue
+
+        # Parse primal variable values (x variables only)
+        if section == "primal_columns":
+            parts = line.split()
+            if len(parts) == 2:
+                var_id, value = parts[0], parts[1]
+                if (
+                    var_id.startswith("x")
+                    and var_id in symbol_map.bySymbol
+                    and var_id not in ["ONE_VAR_CONSTANT", "x2"]
+                ):
+                    symbol_map.bySymbol[var_id].value = float(value)
+
+        # Parse constraint dual values (c_ constraints only) from dual rows
+        elif section == "dual_rows":
+            parts = line.split()
+            if len(parts) == 2:
+                constraint_id, dual = parts[0], parts[1]
+                if (
+                    constraint_id.startswith("c_")
+                    and constraint_id in symbol_map.bySymbol
+                    and constraint_id != "c_e_ONE_VAR_CONSTANT"
+                ):
+                    instance.dual[symbol_map.bySymbol[constraint_id]] = float(dual)
+
+    results = Results(
+        solver_status=solver_status, termination_condition=termination_condition
+    )
+
+    return instance, results, dynamic_components
+
+
 def load_problem_info(prob_sol_files_directory):
     with open(
         os.path.join(prob_sol_files_directory, "instance.pickle"), "rb"