Skip to content

Commit

Permalink
Merge branch 'main' into dispatches_cleanup
Browse files Browse the repository at this point in the history
  • Loading branch information
adam-a-a authored Dec 21, 2023
2 parents 0d93d50 + a9e9b7a commit 0b52a0e
Show file tree
Hide file tree
Showing 53 changed files with 38,373 additions and 45,483 deletions.
8 changes: 7 additions & 1 deletion .readthedocs.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,13 @@ version: 2
sphinx:
fail_on_warning: true

build:
os: "ubuntu-22.04"
tools:
python: "3.8"

python:
version: "3.8"
install:
- requirements: requirements-dev.txt
- method: pip
path: .
Original file line number Diff line number Diff line change
Expand Up @@ -401,7 +401,8 @@
"print(\"Cost Parameters:\")\n",
"print(f\"Wind Capital Cost\\t\\t\\t{wind_cap_cost} $/kW\")\n",
"print(f\"Wind Fixed Operating Cost\\t\\t{wind_op_cost} $/kW-yr\")\n",
"print(f\"Battery Capital Cost\\t\\t\\t{batt_cap_cost} $/kW\")\n",
"print(f\"Battery Power Capital Cost \\t\\t\\t{batt_cap_cost_kw} $/kW\")\n",
"print(f\"Battery Energy Capital Cost \\t\\t\\t{batt_cap_cost_kwh} $/kWh\")\n",
"print(f\"PEM Capital Cost\\t\\t\\t{pem_cap_cost} $/kW\")\n",
"print(f\"PEM Fixed Operating Cost\\t\\t{pem_op_cost} $/kW-yr\")\n",
"print(f\"PEM Variable Operating Cost\\t\\t{pem_var_cost} $/kWh\")\n",
Expand Down Expand Up @@ -829,8 +830,7 @@
"tank_type = 'simple'\n",
"input_params['wind_resource'] = {t:\n",
" {'wind_resource_config': {\n",
" 'resource_speed': [wind_speeds[t]]\n",
" }\n",
" 'capacity_factor': [wind_cfs[t]]}\n",
" } for t in range(7)}\n",
" \n",
"mp_model = MultiPeriodModel(n_time_points=7,\n",
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -1865,7 +1865,6 @@
"source": [
"pd.set_option('display.max_columns', None)\n",
"pd.set_option('display.max_rows', None)\n",
"plt.style.use('seaborn')\n",
"\n",
"df = double_loop_outputs_for_gen(output_dir, rts_gmlc_source_data_path)\n",
"prescient_df = df[df['Model'] == \"Prescient\"].copy()\n",
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -91,7 +91,7 @@ def compute_day_ahead_bids(self, date, hour=0):
return full_bids

def compute_real_time_bids(
self, date, hour, _, __
self, date, hour, realized_day_ahead_prices, realized_day_ahead_dispatches
):
"""
RT Bid: from 0 MW to (Wind Resource - PEM capacity) MW, bid $0/MWh.
Expand Down
26 changes: 26 additions & 0 deletions dispatches/case_studies/renewables_case/README.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,26 @@
# Renewable Energy Case Study: Wind + PEM and Wind + Battery Plants

This directory contains the files required for the renewable energy case studies, which use flowsheets in which any of the following
technologies (modeled as unit models) may be present: Wind, PV, Battery, PEM Electrolysis, Hydrogen Storage Tank and Hydrogen Turbine.
The RE case studies focused on a Wind + PEM plant and a Wind + Battery plant, while the the industry-partnership case study looked at a PV + Battery + PEM + Hydrogen.
Each of these three case studies has different approaches to modeling prices or the grid.

## Wind + PEM case:

1. Price-taker Design Optimization: `dispatches/case_studies/renewables_case/run_pricetaker_wind_PEM.py`
2. Market Surrogate Design Optimization: `dispatches/case_studies/renewables_case/RE_surrogate_optimization_steadystate.py`
3. Double Loop Simulation for Validation: `dispatches/case_studies/renewables_case/run_double_loop_PEM.py`
4. Market Surrogate Design and Validation Plotting: `dispatches/case_studies/renewables_case/SurrogateDesignResults.ipynb`

## Wind + Battery case:
1. Price-taker Design Optimization: `dispatches/case_studies/renewables_case/run_pricetaker_battery_ratio_size.py`
2. Double Loop Simulation: `dispatches/case_studies/renewables_case/run_double_loop_battery.py`
3. Parametrized Bidder Double Loop Simulation: `dispatches/case_studies/renewables_case/run_double_loop_battery_parametrized.py`

## PV + Battery + Hydrogen case:
1. Price-taker Design Optimization: `dispatches/case_studies/renewables_case/solar_battery_hydrogen.py`

## Software and Hardware
All the models are run on a Red Hat Enterprise Linux Server version 7.9 (Maipo). The versions for the solvers used are given below:
- Xpress: Version 8.13.0
- IPOPT: Version 3.13.2
7 changes: 4 additions & 3 deletions dispatches/case_studies/renewables_case/RE_flowsheet.py
Original file line number Diff line number Diff line change
Expand Up @@ -126,8 +126,9 @@ def add_pem(m, outlet_pressure_bar):

m.fs.pem = PEM_Electrolyzer(property_package=m.fs.h2ideal_props)

# Conversion of kW to mol/sec of H2. (elec*elec_to_mol) based on H-tec design of 54.517kW-hr/kg
m.fs.pem.electricity_to_mol.fix(0.002527406)
# Conversion of kW to mol/sec of H2. (elec*elec_to_mol)
# 50 kWh/kg since that is the value for most NEL PEM electrolyzers: https://nelhydrogen.com/product/m-series-3/ (see M3000 data)
m.fs.pem.electricity_to_mol.fix(0.00275984)
m.fs.pem.outlet.pressure.setub(max_pressure_bar * 1e5)
m.fs.pem.outlet.pressure.fix(outlet_pressure_bar * 1e5)
m.fs.pem.outlet.temperature.fix(pem_temp)
Expand All @@ -151,7 +152,7 @@ def add_battery(m, batt_mw):
m.fs.battery.discharging_eta.set_value(0.95)
m.fs.battery.dt.set_value(timestep_hrs)
m.fs.battery.nameplate_power.fix(batt_mw * 1e3)
m.fs.battery.duration = Param(default=4, mutable=True, units=pyunits.kWh/pyunits.kW)
m.fs.battery.duration = Param(default=duration, mutable=True, units=pyunits.kWh/pyunits.kW)
m.fs.battery.four_hr_battery = Constraint(expr=m.fs.battery.nameplate_power * m.fs.battery.duration == m.fs.battery.nameplate_energy)
return m.fs.battery

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,11 @@
initialize_mp, wind_battery_pem_model, wind_battery_pem_mp_block


# RT market only or Both RT and DA markets
rt_market_only = True
include_wind_capital_cost = False
shortfall = 1000

# path for folder that has surrogate models
re_nn_dir = Path(__file__).parent / "data" / "steady_state_surrogate"

Expand All @@ -57,16 +62,21 @@ def load_surrogate_model(re_nn_dir):
pem_clusters_mean = centers[:, 1]
resource_clusters_mean = centers[:, 2]

with open(re_nn_dir / "revenue" / "RE_revenue_params_2_25.json", 'r') as f:
rev_data = json.load(f)

# load keras neural networks
# Input variables are PEM bid price, PEM MW, Reserve Factor and Load Shed Price
nn_rev = keras.models.load_model(re_nn_dir / "revenue" / "RE_revenue_2_25")
nn_dispatch = keras.models.load_model(re_nn_dir / "dispatch_frequency" / "ss_surrogate_model_wind_pmax")

with open(re_nn_dir / "dispatch_frequency" / "ss_surrogate_param_wind_pmax.json", 'r') as f:
dispatch_data = json.load(f)
nn_dispatch = keras.models.load_model(re_nn_dir / "dispatch_frequency" / "ss_surrogate_model_wind_pmax")

if rt_market_only:
rev_data_f = re_nn_dir / "rt_revenue" / "RE_RT_revenue_params_2_25.json"
nn_rev = keras.models.load_model(re_nn_dir / "rt_revenue" / "RE_RT_revenue_2_25")
else:
rev_data_f = re_nn_dir / "revenue" / "RE_revenue_params_2_25.json"
nn_rev = keras.models.load_model(re_nn_dir / "revenue" / "RE_revenue_2_25")

with open(rev_data_f, 'rb') as f:
rev_data = json.load(f)

# load keras models and create OMLT NetworkDefinition objects
#Revenue model definition
Expand Down Expand Up @@ -97,10 +107,10 @@ def conceptual_design_dynamic_RE(input_params, PEM_bid=None, PEM_MW=None, verbos
# add surrogate input to the model
m.wind_system_capacity = Var(domain=NonNegativeReals, bounds=(100 * 1e3, 1000 * 1e3), initialize=input_params['wind_mw'] * 1e3)

m.pem_system_capacity = Var(domain=NonNegativeReals, bounds=(127.05 * 1e3, 423.5 * 1e3), initialize=input_params['pem_mw'] * 1e3, units=pyunits.kW)
m.pem_system_capacity = Var(domain=NonNegativeReals, bounds=(127.5 * 1e3, 423.5 * 1e3), initialize=input_params['pem_mw'] * 1e3, units=pyunits.kW)
m.pem_bid = Var(within=NonNegativeReals, bounds=(15, 45), initialize=45) # Energy Bid $/MWh
m.reserve_percent = Param(within=NonNegativeReals, initialize=15) # Reserves Fraction on Grid
m.shortfall_price = Param(within=NonNegativeReals, initialize=1000) # Energy price during load shed
m.shortfall_price = Param(within=NonNegativeReals, initialize=shortfall) # Energy price during load shed

inputs = [m.pem_bid, m.pem_system_capacity * 1e-3 / 847 * m.wind_system_capacity * 1e-3, m.reserve_percent, m.shortfall_price]

Expand Down Expand Up @@ -213,7 +223,7 @@ def conceptual_design_dynamic_RE(input_params, PEM_bid=None, PEM_MW=None, verbos
scenario_models.append(scenario_model)

m.plant_cap_cost = Expression(
expr=input_params['wind_cap_cost'] * m.wind_system_capacity + input_params['pem_cap_cost'] * m.pem_system_capacity)
expr=input_params['wind_cap_cost'] * m.wind_system_capacity * int(include_wind_capital_cost) + input_params['pem_cap_cost'] * m.pem_system_capacity)
m.annual_fixed_cost = pyo.Expression(
expr=m.wind_system_capacity * input_params["wind_op_cost"] + m.pem_system_capacity * input_params["pem_op_cost"])
m.plant_operation_cost = Expression(
Expand All @@ -222,7 +232,8 @@ def conceptual_design_dynamic_RE(input_params, PEM_bid=None, PEM_MW=None, verbos
expr=sum(scenario_models[i].hydrogen_revenue for i in range(num_rep_days)))

m.NPV = Expression(expr=-m.plant_cap_cost + PA * (m.rev + m.hydrogen_rev - m.plant_operation_cost - m.annual_fixed_cost))
m.obj = Objective(expr=-m.NPV * 1e-8)
m.NPV_ann = Expression(expr=-m.plant_cap_cost / PA + (m.rev + m.hydrogen_rev - m.plant_operation_cost - m.annual_fixed_cost))
m.obj = Objective(expr=-m.NPV * 1e-7)

return m, num_rep_days

Expand Down Expand Up @@ -250,7 +261,8 @@ def record_result(m, num_rep_days):
"pem_bid": value(m.pem_bid),
"e_revenue": value(m.rev),
"h_revenue": value(m.hydrogen_rev),
"NPV": value(m.NPV)
"NPV": value(m.NPV),
"NPV_ann": value(m.NPV_ann)
}

for day in range(num_rep_days):
Expand All @@ -263,6 +275,7 @@ def record_result(m, num_rep_days):
print("Plant Hydrogen Revenue Annual = ${}".format(value(m.hydrogen_rev)))
print("Plant Total Revenue Annual = ${}".format(value(m.rev + m.hydrogen_rev)))
print("Plant NPV = ${}".format(value(m.NPV)))
print("Plant NPV Annualized = ${}".format(value(m.NPV_ann)))

print('----------')
for i in range(num_rep_days):
Expand All @@ -279,8 +292,6 @@ def run_design(PEM_bid=None, PEM_size=None):
model, n_rep_days = conceptual_design_dynamic_RE(default_input_params, PEM_bid=PEM_bid, PEM_MW=PEM_size, verbose=False)
nlp_solver = SolverFactory('ipopt')
nlp_solver.options['max_iter'] = 8000
nlp_solver.options['acceptable_tol'] = 1e-8
nlp_solver.options['bound_push'] = 1e-9
res = nlp_solver.solve(model, tee=True)
if res.Solver.status != 'ok':
solve_log = idaeslog.getInitLogger("infeasibility", idaeslog.INFO, tag="properties")
Expand Down Expand Up @@ -314,7 +325,7 @@ def run_design(PEM_bid=None, PEM_size=None):
"design_opt": True,
"extant_wind": True, # fixed because parameter sweeps didn't change wind size

"wind_cap_cost": wind_cap_cost,
"wind_cap_cost": wind_cap_cost if include_wind_capital_cost else 0,
"wind_op_cost": wind_op_cost,
"pem_cap_cost": pem_cap_cost,
"pem_op_cost": pem_op_cost,
Expand All @@ -323,8 +334,8 @@ def run_design(PEM_bid=None, PEM_size=None):


if __name__ == "__main__":
# result = run_design()
# exit()
result = run_design()
exit()

import multiprocessing as mp
from itertools import product
Expand All @@ -333,8 +344,8 @@ def run_design(PEM_bid=None, PEM_size=None):
sizes = np.linspace(127.05, 423.5, 15)
inputs = product(bids, sizes)

with mp.Pool(processes=4) as p:
with mp.Pool(processes=24) as p:
res = p.starmap(run_design, inputs)

df = pd.DataFrame(res)
df.to_csv("surrogate_results_ss.csv")
df.to_csv(f"wind_PEM/surrogate_results_ss_rt_{shortfall}.csv")
1,593 changes: 933 additions & 660 deletions dispatches/case_studies/renewables_case/SurrogateDesignResults.ipynb

Large diffs are not rendered by default.

124 changes: 124 additions & 0 deletions dispatches/case_studies/renewables_case/battery_parametrized_bidder.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,124 @@
#################################################################################
# DISPATCHES was produced under the DOE Design Integration and Synthesis Platform
# to Advance Tightly Coupled Hybrid Energy Systems program (DISPATCHES), and is
# copyright (c) 2020-2023 by the software owners: The Regents of the University
# of California, through Lawrence Berkeley National Laboratory, National
# Technology & Engineering Solutions of Sandia, LLC, Alliance for Sustainable
# Energy, LLC, Battelle Energy Alliance, LLC, University of Notre Dame du Lac, et
# al. All rights reserved.
#
# Please see the files COPYRIGHT.md and LICENSE.md for full copyright and license
# information, respectively. Both files are also available online at the URL:
# "https://github.com/gmlc-dispatches/dispatches".
#################################################################################
import numpy as np
from idaes.apps.grid_integration.bidder import convert_marginal_costs_to_actual_costs, tx_utils
from dispatches.workflow.parametrized_bidder import ParametrizedBidder


class FixedParametrizedBidder(ParametrizedBidder):

"""
Template class for bidders that use fixed parameters.
The functions for computing the day ahead and real time bids do not use any information from Prescient,
and only depend on internal system information such as wind capacity factors, and on bid parameters.
"""

def __init__(
self,
bidding_model_object,
day_ahead_horizon,
real_time_horizon,
solver,
forecaster,
storage_marginal_cost,
storage_mw
):
super().__init__(bidding_model_object,
day_ahead_horizon,
real_time_horizon,
solver,
forecaster)
self.wind_marginal_cost = 0
self.wind_mw = self.bidding_model_object._wind_pmax_mw
self.storage_marginal_cost = storage_marginal_cost
self.storage_mw = storage_mw

def compute_day_ahead_bids(self, date, hour=0):
"""
Day ahead bid has two parts:
1. the part of the DA wind energy that is able to be stored in the battery is at the higher "storage_marginal_cost"
2. the remainder of the wind energy is bid at the wind marginal cost of $0/MWh
For each time period in the day ahead horizon, the marginal cost bid is assembled as the two parts.
Then the marginal costs are converted to actual costs, as expected by Prescient.
"""
gen = self.generator
forecast = self.forecaster.forecast_day_ahead_capacity_factor(date, hour, gen, self.day_ahead_horizon)

full_bids = {}

for t_idx in range(self.day_ahead_horizon):
da_wind = forecast[t_idx] * self.wind_mw
p_max = max(da_wind, self.storage_mw)
bids = [(0, 0), (max(0, da_wind - self.storage_mw), 0), (p_max, self.storage_marginal_cost)]
cost_curve = convert_marginal_costs_to_actual_costs(bids)

temp_curve = {
"data_type": "cost_curve",
"cost_curve_type": "piecewise",
"values": cost_curve,
}
tx_utils.validate_and_clean_cost_curve(
curve=temp_curve,
curve_type="cost_curve",
p_min=0,
p_max=max([p[0] for p in cost_curve]),
gen_name=gen,
t=t_idx,
)

t = t_idx + hour
full_bids[t] = {}
full_bids[t][gen] = {}
full_bids[t][gen]["p_cost"] = cost_curve
full_bids[t][gen]["p_min"] = 0
full_bids[t][gen]["p_max"] = p_max
full_bids[t][gen]["startup_capacity"] = p_max
full_bids[t][gen]["shutdown_capacity"] = p_max

self._record_bids(full_bids, date, hour, Market="Day-ahead")
return full_bids

def compute_real_time_bids(
self, date, hour, realized_day_ahead_prices, realized_day_ahead_dispatches
):
"""
Real time bid has two parts:
1. the part of the RT wind energy that is able to be stored in the battery is at the higher "storage_marginal_cost"
2. the remainder of the wind energy is bid at the wind marginal cost of $0/MWh
For each time period in the day ahead horizon, the marginal cost bid is assembled as the two parts.
Then the marginal costs are converted to actual costs, as expected by Prescient.
"""
gen = self.generator
forecast = self.forecaster.forecast_real_time_capacity_factor(date, hour, gen, self.day_ahead_horizon)

full_bids = {}

for t_idx in range(self.real_time_horizon):
rt_wind = forecast[t_idx] * self.wind_mw
p_max = max(rt_wind, self.storage_mw)
bids = [(0, 0), (max(0, rt_wind - self.storage_mw), 0), (p_max, self.storage_marginal_cost)]

t = t_idx + hour
full_bids[t] = {}
full_bids[t][gen] = {}
full_bids[t][gen]["p_cost"] = convert_marginal_costs_to_actual_costs(bids)
full_bids[t][gen]["p_min"] = 0
full_bids[t][gen]["p_max"] = p_max
full_bids[t][gen]["startup_capacity"] = p_max
full_bids[t][gen]["shutdown_capacity"] = p_max

self._record_bids(full_bids, date, hour, Market="Real-time")
return full_bids
Binary file not shown.
Binary file not shown.
Loading

0 comments on commit 0b52a0e

Please sign in to comment.