diff --git a/examples/11-load-following-long-duration-energy-storage.ipynb b/examples/11-load-following-long-duration-energy-storage.ipynb index 3cc425e0e..b72d64fd3 100644 --- a/examples/11-load-following-long-duration-energy-storage.ipynb +++ b/examples/11-load-following-long-duration-energy-storage.ipynb @@ -21,7 +21,7 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "metadata": {}, "outputs": [ { diff --git a/examples/inputs/11-load-following-long-duration-energy-storage.yaml b/examples/inputs/11-load-following-long-duration-energy-storage.yaml index af6102e5b..bd092183e 100644 --- a/examples/inputs/11-load-following-long-duration-energy-storage.yaml +++ b/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 diff --git a/examples/legacy/CSP_PV_Battery_Analysis/simulation_init.py b/examples/legacy/CSP_PV_Battery_Analysis/simulation_init.py index 613973861..fd280b3da 100644 --- a/examples/legacy/CSP_PV_Battery_Analysis/simulation_init.py +++ b/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 }, diff --git a/hopp/simulation/technologies/dispatch/hybrid_dispatch_builder_solver.py b/hopp/simulation/technologies/dispatch/hybrid_dispatch_builder_solver.py index 725297547..ea6cb4d64 100644 --- a/hopp/simulation/technologies/dispatch/hybrid_dispatch_builder_solver.py +++ b/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( diff --git a/hopp/simulation/technologies/dispatch/hybrid_dispatch_options.py b/hopp/simulation/technologies/dispatch/hybrid_dispatch_options.py index d19a78e72..86165db22 100644 --- a/hopp/simulation/technologies/dispatch/hybrid_dispatch_options.py +++ b/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 diff --git a/pyproject.toml b/pyproject.toml index 02a1522fc..77aac9905 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -20,6 +20,7 @@ dependencies = [ "future", "global_land_mask", "hybridbosse", + "highspy", "lcoe", "lxml", "matplotlib", diff --git a/tests/hopp/test_battery_dispatch.py b/tests/hopp/test_battery_dispatch.py index a83ddc046..ca278a95b 100644 --- a/tests/hopp/test_battery_dispatch.py +++ b/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"): diff --git a/tests/hopp/test_dispatch.py b/tests/hopp/test_dispatch.py index 27fd715e7..f8ed2696b 100644 --- a/tests/hopp/test_dispatch.py +++ b/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 diff --git a/tests/hopp/test_generic_plant.py b/tests/hopp/test_generic_plant.py index 0b489c929..6f1bf9793 100644 --- a/tests/hopp/test_generic_plant.py +++ b/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,