Update default solver

examples/11-load-following-long-duration-energy-storage.ipynb

@@ -21,7 +21,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": null,
    "metadata": {},
    "outputs": [
     {

examples/inputs/11-load-following-long-duration-energy-storage.yaml

@@ -36,7 +36,7 @@ technologies:
     turbine_rating_kw: 5000
     fin_model: !include default_fin_config.yaml
   battery: # VRDB
-    system_capacity_kwh: 100000
+    system_capacity_kwh: 200000
     system_capacity_kw: 10000
     minimum_SOC: 20.0
     maximum_SOC: 100.0

examples/legacy/CSP_PV_Battery_Analysis/simulation_init.py

@@ -142,7 +142,7 @@ def init_hybrid_plant(techs_in_sim: list, is_test: bool = False, ud_techs: dict
                                     dispatch_options={
                                         'is_test_start_year': is_test,
                                         'is_test_end_year': is_test,
-                                        'solver': 'cbc',
+                                        'solver': 'appsi_highs',
                                         'grid_charging': False,
                                         'pv_charging_only': True
                                         },

hopp/simulation/technologies/dispatch/hybrid_dispatch_builder_solver.py

@@ -139,6 +139,8 @@ def solve_dispatch_model(self, start_time: int, n_days: int):
             solver_results = self.glpk_solve()
         elif self.options.solver == "cbc":
             solver_results = self.cbc_solve()
+        elif self.options.solver == "highs":
+            solver_results = self.highs_solve()
         elif self.options.solver == "xpress":
             solver_results = self.xpress_solve()
         elif self.options.solver == "xpress_persistent":
@@ -187,6 +189,96 @@ def glpk_solve(self):
         return HybridDispatchBuilderSolver.glpk_solve_call(
             self.pyomo_model, self.options.log_name, self.options.solver_options
         )
+
+    @staticmethod
+    def highs_solve_call(
+        pyomo_model: pyomo.ConcreteModel,
+        log_name: str = "",
+        user_solver_options: dict = None,
+    ):
+
+        # log_name = "annual_solve_GLPK.log"  # For debugging MILP solver
+        # Ref. on solver options: https://en.wikibooks.org/wiki/GLPK/Using_GLPSOL
+        # highs_solver_options = {
+        #     "cuts": None,
+        #     "presol": None,
+        #     # 'mostf': None,
+        #     # 'mipgap': 0.001,
+        #     "tmlim": 30,
+        # }
+        highs_solver_options = dict(
+                time_limit=60.0,
+                mip_rel_gap=0.05,  # TODO ???
+                presolve = 'on',
+                # simplex_strategy = "dual",
+                solver="ipm"
+            )
+
+        solver_options = SolverOptions(
+            highs_solver_options, log_name, user_solver_options, "log"
+        )
+        with pyomo.SolverFactory("appsi_highs") as solver:
+            results = solver.solve(pyomo_model, options=solver_options.constructed)
+        HybridDispatchBuilderSolver.log_and_solution_check(
+            log_name,
+            solver_options.instance_log,
+            results.solver.termination_condition,
+            pyomo_model,
+        )
+        return results
+
+    def highs_solve(self):
+        return HybridDispatchBuilderSolver.highs_solve_call(
+            self.pyomo_model, self.options.log_name, self.options.solver_options
+        )
+
+    @staticmethod
+    def scip_solve_call(
+        pyomo_model: pyomo.ConcreteModel,
+        log_name: str = "",
+        user_solver_options: dict = None,
+    ):
+
+        # log_name = "annual_solve_GLPK.log"  # For debugging MILP solver
+        # Ref. on solver options: https://en.wikibooks.org/wiki/GLPK/Using_GLPSOL
+        # highs_solver_options = {
+        #     "cuts": None,
+        #     "presol": None,
+        #     # 'mostf': None,
+        #     # 'mipgap': 0.001,
+        #     "tmlim": 30,
+        # }
+        scip_solver_options = {
+                "limits/gap": 0.5,
+                "limits/time": 60.0,
+                "display/freq": 0.5,
+                "presolving/maxrounds": -1,
+                "limits/nodes": 1000000,
+                "heuristics/emphasis": "aggressive",
+                "separating/maxrounds": 10,
+                "separating/maxcuts": 100,
+                "branching/priority": "fullstrong"
+                # this is currently useless, as pyomo is not calling the concurrent solver
+                # 'parallel/maxnthreads': 16,
+            }
+
+        solver_options = SolverOptions(
+            scip_solver_options, log_name, user_solver_options, "log"
+        )
+        with pyomo.SolverFactory("scip") as solver:
+            results = solver.solve(pyomo_model, options=solver_options.constructed)
+        HybridDispatchBuilderSolver.log_and_solution_check(
+            log_name,
+            solver_options.instance_log,
+            results.solver.termination_condition,
+            pyomo_model,
+        )
+        return results
+
+    def scip_solve(self):
+        return HybridDispatchBuilderSolver.scip_solve_call(
+            self.pyomo_model, self.options.log_name, self.options.solver_options
+        )
 
     @staticmethod
     def gurobi_ampl_solve_call(

hopp/simulation/technologies/dispatch/hybrid_dispatch_options.py

@@ -17,7 +17,7 @@ class HybridDispatchOptions:
     Args:
         dispatch_options (dict): Contains attribute key-value pairs to change default options.
 
-            - **solver** (str, default='cbc'): MILP solver used for dispatch optimization problem. Options are `('glpk', 'cbc', 'xpress', 'xpress_persistent', 'gurobi_ampl', 'gurobi')`.
+            - **solver** (str, default='highs'): MILP solver used for dispatch optimization problem. Options are `('glpk', 'cbc', 'highs', 'scip', 'xpress', 'xpress_persistent', 'gurobi_ampl', 'gurobi')`.
 
             - **solver_options** (dict): Dispatch solver options.
 
@@ -60,7 +60,7 @@ class HybridDispatchOptions:
     """
 
     def __init__(self, dispatch_options: dict = None):
-        self.solver: str = "cbc"
+        self.solver: str = "highs"
         self.solver_options: dict = (
             {}
         )  # used to update solver options, look at specific solver for option names

pyproject.toml

@@ -20,6 +20,7 @@ dependencies = [
     "future",
     "global_land_mask",
     "hybridbosse",
+    "highspy",
     "lcoe",
     "lxml",
     "matplotlib",

tests/hopp/test_battery_dispatch.py

@@ -15,6 +15,8 @@
 from hopp import ROOT_DIR
 from tests.hopp.utils import DEFAULT_FIN_CONFIG
 
+import matplotlib.pyplot as plt
+
 solar_resource_file = ROOT_DIR / "simulation" / "resource_files" / "solar" / "35.2018863_-101.945027_psmv3_60_2012.csv"
 wind_resource_file = ROOT_DIR / "simulation" / "resource_files" / "wind" / "35.2018863_-101.945027_windtoolkit_2012_60min_80m_100m.srw"
 site = SiteInfo(flatirons_site, solar_resource_file=solar_resource_file, wind_resource_file=wind_resource_file)
@@ -57,36 +59,90 @@ def create_test_objective_rule(m):
 
 def test_batterystateless_dispatch(subtests):
     expected_objective = 28957.15
-
-    # Run battery stateful as system model first
-    technologies = technologies_input.copy()
-    technologies['battery']['tracking'] = True
-    model = pyomo.ConcreteModel(name='battery_only')
-    model.forecast_horizon = pyomo.Set(initialize=range(dispatch_n_look_ahead))
-    model.price = pyomo.Param(model.forecast_horizon,
-                              within=pyomo.Reals,
-                              initialize=prices,
-                              mutable=True,
-                              units=u.USD / u.MWh)
-
-    config = BatteryConfig.from_dict(technologies['battery'])
-    battery = Battery(site, config=config)
-    battery._dispatch = SimpleBatteryDispatch(model,
-                                              model.forecast_horizon,
-                                              battery._system_model,
-                                              battery._financial_model,
-                                              'battery',
-                                              HybridDispatchOptions())
+    solvers = {#'glpk': HybridDispatchBuilderSolver.glpk_solve_call, 
+               'highs': HybridDispatchBuilderSolver.highs_solve_call, 
+               #'cbc': HybridDispatchBuilderSolver.cbc_solve_call, 
+               #'scip': HybridDispatchBuilderSolver.scip_solve_call
+               }
+    import time, copy
+    fig, ax = plt.subplots(len(solvers), 1, figsize=(10, 5), sharex=True)
+    fig_single, ax_single = plt.subplots(1, 1, figsize=(10, 5))
+    line_types = ['-', '--', '-.', ':']
+    print(f"Expected objective: {expected_objective:.2f}", flush=True)
+    for i, (solver_name, results) in enumerate(solvers.items()):
+        # Run battery stateful as system model first
+        technologies = technologies_input.copy()
+        technologies['battery']['tracking'] = True
+        model = pyomo.ConcreteModel(name='battery_only')
+        model.forecast_horizon = pyomo.Set(initialize=range(dispatch_n_look_ahead))
+        model.price = pyomo.Param(model.forecast_horizon,
+                                within=pyomo.Reals,
+                                initialize=prices,
+                                mutable=True,
+                                units=u.USD / u.MWh)
+
+        config = BatteryConfig.from_dict(technologies['battery'])
+        battery = Battery(site, config=config)
+        solver = "notsolver"
+        battery._dispatch = SimpleBatteryDispatch(model,
+                                                model.forecast_horizon,
+                                                battery._system_model,
+                                                battery._financial_model,
+                                                'battery',
+                                                HybridDispatchOptions({"solver": solver}))
+
+        model.test_objective = pyomo.Objective(
+            rule=create_test_objective_rule,
+            sense=pyomo.maximize)
+
+        battery.dispatch.initialize_parameters()
+        battery.dispatch.update_time_series_parameters(0)
+        battery.dispatch.update_dispatch_initial_soc(battery.dispatch.minimum_soc)   # Set initial SOC to minimum
+        assert_units_consistent(model)
 
-    model.test_objective = pyomo.Objective(
-        rule=create_test_objective_rule,
-        sense=pyomo.maximize)
-
-    battery.dispatch.initialize_parameters()
-    battery.dispatch.update_time_series_parameters(0)
-    battery.dispatch.update_dispatch_initial_soc(battery.dispatch.minimum_soc)   # Set initial SOC to minimum
-    assert_units_consistent(model)
-    results = HybridDispatchBuilderSolver.glpk_solve_call(model)
+        if False: # for solver debugging
+            t1 = time.time()
+            n = 50
+            battery_dispatch_sum = 0
+            # print(f"Running Solver {i}: {solver_name}", flush=True)
+            for j in range(n):
+                model_internal = model
+
+                results = solvers[solver_name](model_internal)
+
+                battery_dispatch = np.array(battery.dispatch.power)[0:dispatch_n_look_ahead]
+                battery_dispatch_sum += np.sum(battery_dispatch)
+                if j == 0:
+                    if i == 0:
+                        label_dispatch = "dispatch"
+                        label_actual = "actual"
+                    else:
+                        label_dispatch = None
+                        label_actual = None
+                    ax[i].plot(battery.dispatch.power[0:dispatch_n_look_ahead], label=label_dispatch)
+                    ax[i].plot(battery.generation_profile[0:dispatch_n_look_ahead], label=label_actual)
+                    ax[i].set(ylabel='Power (MW)', title=f"Solver: {solver_name}")
+                    ax_single.plot(battery.dispatch.power[0:dispatch_n_look_ahead], linestyle=line_types[i], label=f"{solver_name}")
+                    if i == 0:
+                        ax[i].legend()
+                    if i == len(solvers) - 1:
+                        ax[i].set(xlabel='Time (hours)')
+            t2 = time.time()
+            battery_dispatch_sum /= n
+            print(f"Solver {i}: {solver_name}, N: {n}, TTime: {t2 - t1:.4f} (s), Time/run: {(t2 - t1) / n:.4f} (s), Obj: {pyomo.value(model_internal.test_objective)}, Ave dis. power: {battery_dispatch_sum / n:.4f} (MWh), Exit: {results.solver.termination_condition }", flush=True)
+    if True: # for solver debugging
+        ax_single.set(ylabel='Power (MW)', xlabel='Time (hours)', title="Dispatch Power")
+        ax_single.legend()
+        fig.tight_layout()
+        fig_single.tight_layout()
+        fig.subplots_adjust(hspace=0.5)
+        fig.savefig("solver_comparison.png")
+    fig_single.savefig("solver_comparison_single_plot.png")
+    # assert False
+    # results = HybridDispatchBuilderSolver.glpk_solve_call(model)
+    results = HybridDispatchBuilderSolver.highs_solve_call(model)
+    # results = HybridDispatchBuilderSolver.scip_solve_call(model)
+    # results = HybridDispatchBuilderSolver.cbc_solve_call(model)
 
     with subtests.test("TerminationCondition"):
         assert results.solver.termination_condition == TerminationCondition.optimal
@@ -144,7 +200,8 @@ def test_batterystateless_dispatch(subtests):
     battery_sl.dispatch.initialize_parameters()
     battery_sl.dispatch.update_time_series_parameters(0)
     assert_units_consistent(model_sl)
-    results = HybridDispatchBuilderSolver.glpk_solve_call(model_sl)
+    # results = HybridDispatchBuilderSolver.glpk_solve_call(model_sl)
+    results = HybridDispatchBuilderSolver.highs_solve_call(model_sl)
 
     with subtests.test("sum_charge_power"):
         assert results.solver.termination_condition == TerminationCondition.optimal
@@ -164,22 +221,22 @@ def test_batterystateless_dispatch(subtests):
             dispatch_power = battery_sl.dispatch.power[i] * 1e3
             assert battery_sl.outputs.P[i] == pytest.approx(dispatch_power, 1e-3 * abs(dispatch_power))
 
-    battery_dispatch = np.array(battery.dispatch.power)[0:48]
+    battery_dispatch = np.array(battery.dispatch.power)[0:dispatch_n_look_ahead]
     battery_actual = np.array(battery.generation_profile[0:dispatch_n_look_ahead]) * 1e-3   # convert to MWh
-    battery_sl_dispatch = np.array(battery_sl.dispatch.power)[0:48]
-    battery_sl_actual = np.array(battery_sl.generation_profile)[0:48] * 1e-3   # convert to MWh
+    battery_sl_dispatch = np.array(battery_sl.dispatch.power)[0:dispatch_n_look_ahead]
+    battery_sl_actual = np.array(battery_sl.generation_profile)[0:dispatch_n_look_ahead] * 1e-3   # convert to MWh
 
     with subtests.test("battery_dispatch vs battery_sl_dispatch"):
-        assert sum(battery_dispatch - battery_sl_dispatch) == 0
+        assert sum(battery_dispatch - battery_sl_dispatch) == pytest.approx(0.0)
 
     with subtests.test("battery_actual vs battery_dispatch"):
-        assert sum(abs(battery_actual - battery_dispatch)) <= 33.5
+        assert sum(abs(battery_actual - battery_dispatch)) <= 33.9
 
     with subtests.test("battery_sl_actual vs battery_sl_dispatch"):
-        assert sum(abs(battery_sl_actual - battery_sl_dispatch)) == 0
+        assert sum(abs(battery_sl_actual - battery_sl_dispatch)) == pytest.approx(0.0)
 
     with subtests.test("battery_actual vs battery_sl_actual"):
-        assert sum(abs(battery_actual - battery_sl_actual)) <= 33.5
+        assert sum(abs(battery_actual - battery_sl_actual)) <= 33.9
 
     with subtests.test("lifecycles_per_day"):
         assert battery_sl.outputs.lifecycles_per_day[0:2] == pytest.approx([0.75048, 1.50096], rel=1e-3)
@@ -214,7 +271,8 @@ def test_batterystateless_cycle_limits(subtests):
     battery_sl.dispatch.initialize_parameters()
     battery_sl.dispatch.update_time_series_parameters(0)
     assert_units_consistent(model_sl)
-    results = HybridDispatchBuilderSolver.glpk_solve_call(model_sl)
+    # results = HybridDispatchBuilderSolver.glpk_solve_call(model_sl)
+    results = HybridDispatchBuilderSolver.highs_solve_call(model_sl)
 
 
     with subtests.test("termination_condition"):

tests/hopp/test_dispatch.py

@@ -105,7 +105,8 @@ def create_test_objective_rule(m):
 
     solar._dispatch.update_time_series_parameters(0)
 
-    results = HybridDispatchBuilderSolver.glpk_solve_call(model)
+    # results = HybridDispatchBuilderSolver.glpk_solve_call(model)
+    results = HybridDispatchBuilderSolver.highs_solve_call(model)
     # results = HybridDispatchBuilderSolver.cbc_solve_call(model)
     # results = HybridDispatchBuilderSolver.xpress_solve_call(model)
     assert results.solver.termination_condition == TerminationCondition.optimal

tests/hopp/test_generic_plant.py

@@ -97,7 +97,7 @@ def hybrid_tech_config():
 def dispatch_options():
     dispatch_opt = {
         "battery_dispatch": "load_following_heuristic",
-        "solver": "cbc",
+        "solver": "highs",
         "n_look_ahead_periods": 48,
         "grid_charging": False,
         "pv_charging_only": False,