From 31437789121ba03c453ea76b24e8d400e5f35303 Mon Sep 17 00:00:00 2001 From: Bo Yuan Date: Wed, 31 Jan 2024 15:27:58 -0500 Subject: [PATCH] get a new set of examples --- ...{ex_opf_wo_renew.py => 01_opf_baseline.py} | 21 +- examples/02_opf_w_esr.py | 133 ++++++++ examples/03_opf_w_cpny.py | 133 ++++++++ examples/04_opf_w_esr_cpny.py | 133 ++++++++ examples/05_opf_w_vre.py | 153 +++++++++ examples/06_opf_w_vre_esr.py | 153 +++++++++ examples/07_opf_w_vre_cpny.py | 153 +++++++++ examples/08_opf_w_vre_esr_cpny.py | 153 +++++++++ examples/09_opf_w_elec_vre.py | 169 ++++++++++ examples/10_opf_w_elec_vre_esr.py | 169 ++++++++++ examples/11_opf_w_elec_vre_cpny.py | 169 ++++++++++ examples/12_opf_w_elec_vre_esr_cpny.py | 169 ++++++++++ examples/ex_opf_w_future_solar.py | 290 ---------------- examples/ex_opf_w_offshore_wind.py | 260 --------------- examples/ex_opf_w_renew.py | 260 --------------- examples/ex_opf_w_renew_elec.py | 310 ------------------ nygrid/run_nygrid.py | 169 +++++++++- 17 files changed, 1867 insertions(+), 1130 deletions(-) rename examples/{ex_opf_wo_renew.py => 01_opf_baseline.py} (85%) create mode 100644 examples/02_opf_w_esr.py create mode 100644 examples/03_opf_w_cpny.py create mode 100644 examples/04_opf_w_esr_cpny.py create mode 100644 examples/05_opf_w_vre.py create mode 100644 examples/06_opf_w_vre_esr.py create mode 100644 examples/07_opf_w_vre_cpny.py create mode 100644 examples/08_opf_w_vre_esr_cpny.py create mode 100644 examples/09_opf_w_elec_vre.py create mode 100644 examples/10_opf_w_elec_vre_esr.py create mode 100644 examples/11_opf_w_elec_vre_cpny.py create mode 100644 examples/12_opf_w_elec_vre_esr_cpny.py delete mode 100644 examples/ex_opf_w_future_solar.py delete mode 100644 examples/ex_opf_w_offshore_wind.py delete mode 100644 examples/ex_opf_w_renew.py delete mode 100644 examples/ex_opf_w_renew_elec.py diff --git a/examples/ex_opf_wo_renew.py b/examples/01_opf_baseline.py similarity index 85% rename from examples/ex_opf_wo_renew.py rename to examples/01_opf_baseline.py index e0026d7..0c2a096 100644 --- a/examples/ex_opf_wo_renew.py +++ b/examples/01_opf_baseline.py @@ -1,7 +1,10 @@ """ -Run multi-period OPF with 2018 data -without renewable generators - +Run multi-period OPF with 2018 data. +1. Baseline case: +No CPNY and CHPE HVDC lines. +No ESRs. +No future solar and offshore wind. +No residential building electrification. """ # %% Packages @@ -19,10 +22,12 @@ # %% Simulation settings # NOTE: Change the following settings to run the simulation - sim_name = 'wo_renew' + sim_name = 'baseline' leading_hours = 12 w_cpny = False # True: add CPNY and CHPE HVDC lines; False: no CPNY and CHPE HVDC lines w_esr = False # True: add ESRs; False: no ESRs + w_vre = False # True: add future solar and offshore wind; False: no future solar and offshore wind + w_elec = False # True: add residential building electrification; False: no residential building electrification start_date = datetime(2018, 1, 1, 0, 0, 0) end_date = datetime(2018, 12, 31, 0, 0, 0) @@ -72,9 +77,9 @@ if w_cpny: grid_data['dcline_prop'] = dcline_prop - logging.info('With CPNY and CHPE HVDC lines.') + logging.info('With CPNY and CHPE HVDC lines.') # Existing and planned HVDC lines else: - grid_data['dcline_prop'] = dcline_prop[4:] + grid_data['dcline_prop'] = dcline_prop[:4] # Only existing HVDC lines logging.info('Without CPNY and CHPE HVDC lines.') # Read ESR property file @@ -83,10 +88,10 @@ if w_esr: logging.info('With ESRs.') - grid_data['esr_prop'] = esr_prop + grid_data['esr_prop'] = esr_prop # Existing and planned ESRs else: logging.info('No ESRs.') - grid_data['esr_prop'] = None + grid_data['esr_prop'] = esr_prop[:8] # Only existing ESRs # %% Set up OPF model diff --git a/examples/02_opf_w_esr.py b/examples/02_opf_w_esr.py new file mode 100644 index 0000000..96e303c --- /dev/null +++ b/examples/02_opf_w_esr.py @@ -0,0 +1,133 @@ +""" +Run multi-period OPF with 2018 data. +2. With ESR case: +No CPNY and CHPE HVDC lines. +Add 6 GW ESRs. +No future solar and offshore wind. +No residential building electrification. +""" +# %% Packages + +import logging +import os +import pickle +import time +from datetime import datetime, timedelta + +import pandas as pd + +from nygrid.run_nygrid import read_grid_data, run_nygrid_one_day + +if __name__ == '__main__': + + # %% Simulation settings + # NOTE: Change the following settings to run the simulation + sim_name = 'w_esr' + leading_hours = 12 + w_cpny = False # True: add CPNY and CHPE HVDC lines; False: no CPNY and CHPE HVDC lines + w_esr = True # True: add ESRs; False: no ESRs + w_vre = False # True: add future solar and offshore wind; False: no future solar and offshore wind + w_elec = False # True: add residential building electrification; False: no residential building electrification + + start_date = datetime(2018, 1, 1, 0, 0, 0) + end_date = datetime(2018, 12, 31, 0, 0, 0) + timestamp_list = pd.date_range(start_date, end_date, freq='1D') + + # %% Set up logging + logging.basicConfig(level=logging.INFO, + format='%(asctime)s - %(levelname)s - %(message)s', + handlers=[logging.FileHandler(f'ex_opf_{sim_name}.log'), + logging.StreamHandler()]) + + prog_start = time.time() + logging.info(f'Start running multi-period OPF simulation {sim_name}.') + + # %% Set up directories + cwd = os.getcwd() + if 'examples' in cwd: + parent_dir = os.path.dirname(cwd) + data_dir = os.path.join(parent_dir, 'data') + else: + data_dir = os.path.join(cwd, 'data') + + grid_data_dir = os.path.join(data_dir, 'grid') + if not os.path.exists(grid_data_dir): + raise FileNotFoundError('Grid data directory not found.') + + logging.info('Grid data directory: {}'.format(grid_data_dir)) + + fig_dir = os.path.join(os.path.dirname(data_dir), 'figures') + logging.info('Figure directory: {}'.format(fig_dir)) + + results_dir = os.path.join(os.path.dirname(data_dir), 'results') + logging.info('Results directory: {}'.format(results_dir)) + + sim_results_dir = os.path.join(results_dir, sim_name) + if not os.path.exists(sim_results_dir): + os.mkdir(sim_results_dir) + + # %% Read grid data + + # Read load and generation profiles + grid_data = read_grid_data(grid_data_dir, start_date.year) + + # Read DC line property file + filename = os.path.join(grid_data_dir, 'dcline_prop.csv') + dcline_prop = pd.read_csv(filename, index_col=0) + + if w_cpny: + grid_data['dcline_prop'] = dcline_prop + logging.info('With CPNY and CHPE HVDC lines.') # Existing and planned HVDC lines + else: + grid_data['dcline_prop'] = dcline_prop[:4] # Only existing HVDC lines + logging.info('Without CPNY and CHPE HVDC lines.') + + # Read ESR property file + filename = os.path.join(grid_data_dir, 'esr_prop.csv') + esr_prop = pd.read_csv(filename, index_col=0) + + if w_esr: + logging.info('With ESRs.') + grid_data['esr_prop'] = esr_prop # Existing and planned ESRs + else: + logging.info('No ESRs.') + grid_data['esr_prop'] = esr_prop[:8] # Only existing ESRs + + # %% Set up OPF model + + # Set options + options = {'UsePTDF': True, + 'solver': 'gurobi', + 'PenaltyForLoadShed': 20_000, + 'PenaltyForBranchMwViolation': 10_000, + 'PenaltyForInterfaceMWViolation': 10_000} + + # No initial condition for the first day + last_gen = None + + # Loop through all days + for d in range(len(timestamp_list) - 1): + t = time.time() + + # Run OPF for one day (24 hours) plus leading hours + # The first day is valid, the leading hours are used to dispatch batteries properly + start_datetime = timestamp_list[d] + end_datetime = start_datetime + timedelta(hours=24 + leading_hours) + + nygrid_results = run_nygrid_one_day(start_datetime, end_datetime, grid_data, grid_data_dir, options, last_gen) + + # Set generator initial condition for the next iteration + last_gen = nygrid_results['PG'].loc[start_datetime].to_numpy().squeeze() + + # Save simulation nygrid_results to pickle + filename = f'nygrid_sim_{sim_name}_{start_datetime.strftime("%Y%m%d%H")}.pkl' + with open(os.path.join(sim_results_dir, filename), 'wb') as f: + pickle.dump(nygrid_results, f) + logging.info(f'Saved simulation nygrid_results in {filename}') + + elapsed = time.time() - t + logging.info(f'Finished running for {start_datetime.strftime("%Y-%m-%d")}. Elapsed time: {elapsed:.2f} seconds') + logging.info('-' * 80) + + tot_elapsed = time.time() - prog_start + logging.info(f"Finished multi-period OPF simulation {sim_name}. Total elapsed time: {tot_elapsed:.2f} seconds") diff --git a/examples/03_opf_w_cpny.py b/examples/03_opf_w_cpny.py new file mode 100644 index 0000000..dd2e8ec --- /dev/null +++ b/examples/03_opf_w_cpny.py @@ -0,0 +1,133 @@ +""" +Run multi-period OPF with 2018 data. +3. With CPNY case: +Add CPNY and CHPE HVDC lines. +No ESRs. +No future solar and offshore wind. +No residential building electrification. +""" +# %% Packages + +import logging +import os +import pickle +import time +from datetime import datetime, timedelta + +import pandas as pd + +from nygrid.run_nygrid import read_grid_data, run_nygrid_one_day + +if __name__ == '__main__': + + # %% Simulation settings + # NOTE: Change the following settings to run the simulation + sim_name = 'w_cpny' + leading_hours = 12 + w_cpny = True # True: add CPNY and CHPE HVDC lines; False: no CPNY and CHPE HVDC lines + w_esr = False # True: add ESRs; False: no ESRs + w_vre = False # True: add future solar and offshore wind; False: no future solar and offshore wind + w_elec = False # True: add residential building electrification; False: no residential building electrification + + start_date = datetime(2018, 1, 1, 0, 0, 0) + end_date = datetime(2018, 12, 31, 0, 0, 0) + timestamp_list = pd.date_range(start_date, end_date, freq='1D') + + # %% Set up logging + logging.basicConfig(level=logging.INFO, + format='%(asctime)s - %(levelname)s - %(message)s', + handlers=[logging.FileHandler(f'ex_opf_{sim_name}.log'), + logging.StreamHandler()]) + + prog_start = time.time() + logging.info(f'Start running multi-period OPF simulation {sim_name}.') + + # %% Set up directories + cwd = os.getcwd() + if 'examples' in cwd: + parent_dir = os.path.dirname(cwd) + data_dir = os.path.join(parent_dir, 'data') + else: + data_dir = os.path.join(cwd, 'data') + + grid_data_dir = os.path.join(data_dir, 'grid') + if not os.path.exists(grid_data_dir): + raise FileNotFoundError('Grid data directory not found.') + + logging.info('Grid data directory: {}'.format(grid_data_dir)) + + fig_dir = os.path.join(os.path.dirname(data_dir), 'figures') + logging.info('Figure directory: {}'.format(fig_dir)) + + results_dir = os.path.join(os.path.dirname(data_dir), 'results') + logging.info('Results directory: {}'.format(results_dir)) + + sim_results_dir = os.path.join(results_dir, sim_name) + if not os.path.exists(sim_results_dir): + os.mkdir(sim_results_dir) + + # %% Read grid data + + # Read load and generation profiles + grid_data = read_grid_data(grid_data_dir, start_date.year) + + # Read DC line property file + filename = os.path.join(grid_data_dir, 'dcline_prop.csv') + dcline_prop = pd.read_csv(filename, index_col=0) + + if w_cpny: + grid_data['dcline_prop'] = dcline_prop + logging.info('With CPNY and CHPE HVDC lines.') # Existing and planned HVDC lines + else: + grid_data['dcline_prop'] = dcline_prop[:4] # Only existing HVDC lines + logging.info('Without CPNY and CHPE HVDC lines.') + + # Read ESR property file + filename = os.path.join(grid_data_dir, 'esr_prop.csv') + esr_prop = pd.read_csv(filename, index_col=0) + + if w_esr: + logging.info('With ESRs.') + grid_data['esr_prop'] = esr_prop # Existing and planned ESRs + else: + logging.info('No ESRs.') + grid_data['esr_prop'] = esr_prop[:8] # Only existing ESRs + + # %% Set up OPF model + + # Set options + options = {'UsePTDF': True, + 'solver': 'gurobi', + 'PenaltyForLoadShed': 20_000, + 'PenaltyForBranchMwViolation': 10_000, + 'PenaltyForInterfaceMWViolation': 10_000} + + # No initial condition for the first day + last_gen = None + + # Loop through all days + for d in range(len(timestamp_list) - 1): + t = time.time() + + # Run OPF for one day (24 hours) plus leading hours + # The first day is valid, the leading hours are used to dispatch batteries properly + start_datetime = timestamp_list[d] + end_datetime = start_datetime + timedelta(hours=24 + leading_hours) + + nygrid_results = run_nygrid_one_day(start_datetime, end_datetime, grid_data, grid_data_dir, options, last_gen) + + # Set generator initial condition for the next iteration + last_gen = nygrid_results['PG'].loc[start_datetime].to_numpy().squeeze() + + # Save simulation nygrid_results to pickle + filename = f'nygrid_sim_{sim_name}_{start_datetime.strftime("%Y%m%d%H")}.pkl' + with open(os.path.join(sim_results_dir, filename), 'wb') as f: + pickle.dump(nygrid_results, f) + logging.info(f'Saved simulation nygrid_results in {filename}') + + elapsed = time.time() - t + logging.info(f'Finished running for {start_datetime.strftime("%Y-%m-%d")}. Elapsed time: {elapsed:.2f} seconds') + logging.info('-' * 80) + + tot_elapsed = time.time() - prog_start + logging.info(f"Finished multi-period OPF simulation {sim_name}. Total elapsed time: {tot_elapsed:.2f} seconds") diff --git a/examples/04_opf_w_esr_cpny.py b/examples/04_opf_w_esr_cpny.py new file mode 100644 index 0000000..d8e60dd --- /dev/null +++ b/examples/04_opf_w_esr_cpny.py @@ -0,0 +1,133 @@ +""" +Run multi-period OPF with 2018 data. +4. With ESR and CPNY case: +Add CPNY and CHPE HVDC lines. +Add 6 GW ESRs. +No future solar and offshore wind. +No residential building electrification. +""" +# %% Packages + +import logging +import os +import pickle +import time +from datetime import datetime, timedelta + +import pandas as pd + +from nygrid.run_nygrid import read_grid_data, run_nygrid_one_day + +if __name__ == '__main__': + + # %% Simulation settings + # NOTE: Change the following settings to run the simulation + sim_name = 'w_esr_cpny' + leading_hours = 12 + w_cpny = True # True: add CPNY and CHPE HVDC lines; False: no CPNY and CHPE HVDC lines + w_esr = True # True: add ESRs; False: no ESRs + w_vre = False # True: add future solar and offshore wind; False: no future solar and offshore wind + w_elec = False # True: add residential building electrification; False: no residential building electrification + + start_date = datetime(2018, 1, 1, 0, 0, 0) + end_date = datetime(2018, 12, 31, 0, 0, 0) + timestamp_list = pd.date_range(start_date, end_date, freq='1D') + + # %% Set up logging + logging.basicConfig(level=logging.INFO, + format='%(asctime)s - %(levelname)s - %(message)s', + handlers=[logging.FileHandler(f'ex_opf_{sim_name}.log'), + logging.StreamHandler()]) + + prog_start = time.time() + logging.info(f'Start running multi-period OPF simulation {sim_name}.') + + # %% Set up directories + cwd = os.getcwd() + if 'examples' in cwd: + parent_dir = os.path.dirname(cwd) + data_dir = os.path.join(parent_dir, 'data') + else: + data_dir = os.path.join(cwd, 'data') + + grid_data_dir = os.path.join(data_dir, 'grid') + if not os.path.exists(grid_data_dir): + raise FileNotFoundError('Grid data directory not found.') + + logging.info('Grid data directory: {}'.format(grid_data_dir)) + + fig_dir = os.path.join(os.path.dirname(data_dir), 'figures') + logging.info('Figure directory: {}'.format(fig_dir)) + + results_dir = os.path.join(os.path.dirname(data_dir), 'results') + logging.info('Results directory: {}'.format(results_dir)) + + sim_results_dir = os.path.join(results_dir, sim_name) + if not os.path.exists(sim_results_dir): + os.mkdir(sim_results_dir) + + # %% Read grid data + + # Read load and generation profiles + grid_data = read_grid_data(grid_data_dir, start_date.year) + + # Read DC line property file + filename = os.path.join(grid_data_dir, 'dcline_prop.csv') + dcline_prop = pd.read_csv(filename, index_col=0) + + if w_cpny: + grid_data['dcline_prop'] = dcline_prop + logging.info('With CPNY and CHPE HVDC lines.') # Existing and planned HVDC lines + else: + grid_data['dcline_prop'] = dcline_prop[:4] # Only existing HVDC lines + logging.info('Without CPNY and CHPE HVDC lines.') + + # Read ESR property file + filename = os.path.join(grid_data_dir, 'esr_prop.csv') + esr_prop = pd.read_csv(filename, index_col=0) + + if w_esr: + logging.info('With ESRs.') + grid_data['esr_prop'] = esr_prop # Existing and planned ESRs + else: + logging.info('No ESRs.') + grid_data['esr_prop'] = esr_prop[:8] # Only existing ESRs + + # %% Set up OPF model + + # Set options + options = {'UsePTDF': True, + 'solver': 'gurobi', + 'PenaltyForLoadShed': 20_000, + 'PenaltyForBranchMwViolation': 10_000, + 'PenaltyForInterfaceMWViolation': 10_000} + + # No initial condition for the first day + last_gen = None + + # Loop through all days + for d in range(len(timestamp_list) - 1): + t = time.time() + + # Run OPF for one day (24 hours) plus leading hours + # The first day is valid, the leading hours are used to dispatch batteries properly + start_datetime = timestamp_list[d] + end_datetime = start_datetime + timedelta(hours=24 + leading_hours) + + nygrid_results = run_nygrid_one_day(start_datetime, end_datetime, grid_data, grid_data_dir, options, last_gen) + + # Set generator initial condition for the next iteration + last_gen = nygrid_results['PG'].loc[start_datetime].to_numpy().squeeze() + + # Save simulation nygrid_results to pickle + filename = f'nygrid_sim_{sim_name}_{start_datetime.strftime("%Y%m%d%H")}.pkl' + with open(os.path.join(sim_results_dir, filename), 'wb') as f: + pickle.dump(nygrid_results, f) + logging.info(f'Saved simulation nygrid_results in {filename}') + + elapsed = time.time() - t + logging.info(f'Finished running for {start_datetime.strftime("%Y-%m-%d")}. Elapsed time: {elapsed:.2f} seconds') + logging.info('-' * 80) + + tot_elapsed = time.time() - prog_start + logging.info(f"Finished multi-period OPF simulation {sim_name}. Total elapsed time: {tot_elapsed:.2f} seconds") diff --git a/examples/05_opf_w_vre.py b/examples/05_opf_w_vre.py new file mode 100644 index 0000000..7f4a740 --- /dev/null +++ b/examples/05_opf_w_vre.py @@ -0,0 +1,153 @@ +""" +Run multi-period OPF with 2018 data. +5. With VRE case: +No CPNY and CHPE HVDC lines. +No ESRs. +Add future solar and offshore wind. +No residential building electrification. +""" +# %% Packages + +import logging +import os +import pickle +import time +from datetime import datetime, timedelta + +import pandas as pd + +from nygrid.run_nygrid import read_grid_data, read_vre_data, run_nygrid_one_day + +if __name__ == '__main__': + + # %% Simulation settings + # NOTE: Change the following settings to run the simulation + sim_name = 'w_vre' + leading_hours = 12 + w_cpny = False # True: add CPNY and CHPE HVDC lines; False: no CPNY and CHPE HVDC lines + w_esr = False # True: add ESRs; False: no ESRs + w_vre = True # True: add future solar and offshore wind; False: no future solar and offshore wind + w_elec = False # True: add residential building electrification; False: no residential building electrification + + start_date = datetime(2018, 1, 1, 0, 0, 0) + end_date = datetime(2018, 12, 31, 0, 0, 0) + timestamp_list = pd.date_range(start_date, end_date, freq='1D') + + # %% Set up logging + logging.basicConfig(level=logging.INFO, + format='%(asctime)s - %(levelname)s - %(message)s', + handlers=[logging.FileHandler(f'ex_opf_{sim_name}.log'), + logging.StreamHandler()]) + + prog_start = time.time() + logging.info(f'Start running multi-period OPF simulation {sim_name}.') + + # %% Set up directories + cwd = os.getcwd() + if 'examples' in cwd: + parent_dir = os.path.dirname(cwd) + data_dir = os.path.join(parent_dir, 'data') + else: + data_dir = os.path.join(cwd, 'data') + + grid_data_dir = os.path.join(data_dir, 'grid') + if not os.path.exists(grid_data_dir): + raise FileNotFoundError('Grid data directory not found.') + + logging.info('Grid data directory: {}'.format(grid_data_dir)) + + fig_dir = os.path.join(os.path.dirname(data_dir), 'figures') + logging.info('Figure directory: {}'.format(fig_dir)) + + results_dir = os.path.join(os.path.dirname(data_dir), 'results') + logging.info('Results directory: {}'.format(results_dir)) + + solar_data_dir = os.path.join(data_dir, 'solar') + print('Solar data directory: {}'.format(solar_data_dir)) + + onshore_wind_data_dir = os.path.join(data_dir, 'onshore_wind') + print('Onshore wind data directory: {}'.format(onshore_wind_data_dir)) + + offshore_wind_data_dir = os.path.join(data_dir, 'offshore_wind') + print('Offshore wind data directory: {}'.format(offshore_wind_data_dir)) + + sim_results_dir = os.path.join(results_dir, sim_name) + if not os.path.exists(sim_results_dir): + os.mkdir(sim_results_dir) + + # %% Read grid data + + # Read load and generation profiles + grid_data = read_grid_data(grid_data_dir, start_date.year) + + # Read DC line property file + filename = os.path.join(grid_data_dir, 'dcline_prop.csv') + dcline_prop = pd.read_csv(filename, index_col=0) + + if w_cpny: + grid_data['dcline_prop'] = dcline_prop + logging.info('With CPNY and CHPE HVDC lines.') # Existing and planned HVDC lines + else: + grid_data['dcline_prop'] = dcline_prop[:4] # Only existing HVDC lines + logging.info('Without CPNY and CHPE HVDC lines.') + + # Read ESR property file + filename = os.path.join(grid_data_dir, 'esr_prop.csv') + esr_prop = pd.read_csv(filename, index_col=0) + + if w_esr: + logging.info('With ESRs.') + grid_data['esr_prop'] = esr_prop # Existing and planned ESRs + else: + logging.info('No ESRs.') + grid_data['esr_prop'] = esr_prop[:8] # Only existing ESRs + + # Read renewable generation profiles + if w_vre: + vre_prop, genmax_profile_vre = read_vre_data(solar_data_dir, onshore_wind_data_dir, offshore_wind_data_dir) + grid_data['vre_prop'] = vre_prop + grid_data['genmax_profile_vre'] = genmax_profile_vre + logging.info('With future solar and offshore wind.') + else: + logging.info('No future solar and offshore wind.') + grid_data['vre_prop'] = None + grid_data['genmax_profile_vre'] = None + + # %% Set up OPF model + + # Set options + options = {'UsePTDF': True, + 'solver': 'gurobi', + 'PenaltyForLoadShed': 20_000, + 'PenaltyForBranchMwViolation': 10_000, + 'PenaltyForInterfaceMWViolation': 10_000} + + # No initial condition for the first day + last_gen = None + + # Loop through all days + for d in range(len(timestamp_list) - 1): + t = time.time() + + # Run OPF for one day (24 hours) plus leading hours + # The first day is valid, the leading hours are used to dispatch batteries properly + start_datetime = timestamp_list[d] + end_datetime = start_datetime + timedelta(hours=24 + leading_hours) + + nygrid_results = run_nygrid_one_day(start_datetime, end_datetime, grid_data, grid_data_dir, options, last_gen) + + # Set generator initial condition for the next iteration + last_gen = nygrid_results['PG'].loc[start_datetime].to_numpy().squeeze() + + # Save simulation nygrid_results to pickle + filename = f'nygrid_sim_{sim_name}_{start_datetime.strftime("%Y%m%d%H")}.pkl' + with open(os.path.join(sim_results_dir, filename), 'wb') as f: + pickle.dump(nygrid_results, f) + logging.info(f'Saved simulation nygrid_results in {filename}') + + elapsed = time.time() - t + logging.info(f'Finished running for {start_datetime.strftime("%Y-%m-%d")}. Elapsed time: {elapsed:.2f} seconds') + logging.info('-' * 80) + + tot_elapsed = time.time() - prog_start + logging.info(f"Finished multi-period OPF simulation {sim_name}. Total elapsed time: {tot_elapsed:.2f} seconds") diff --git a/examples/06_opf_w_vre_esr.py b/examples/06_opf_w_vre_esr.py new file mode 100644 index 0000000..3d0558c --- /dev/null +++ b/examples/06_opf_w_vre_esr.py @@ -0,0 +1,153 @@ +""" +Run multi-period OPF with 2018 data. +6. With ESR and VRE case: +No CPNY and CHPE HVDC lines. +Add 6 GW ESRs. +Add future solar and offshore wind. +No residential building electrification. +""" +# %% Packages + +import logging +import os +import pickle +import time +from datetime import datetime, timedelta + +import pandas as pd + +from nygrid.run_nygrid import read_grid_data, read_vre_data, run_nygrid_one_day + +if __name__ == '__main__': + + # %% Simulation settings + # NOTE: Change the following settings to run the simulation + sim_name = 'w_vre_esr' + leading_hours = 12 + w_cpny = False # True: add CPNY and CHPE HVDC lines; False: no CPNY and CHPE HVDC lines + w_esr = True # True: add ESRs; False: no ESRs + w_vre = True # True: add future solar and offshore wind; False: no future solar and offshore wind + w_elec = False # True: add residential building electrification; False: no residential building electrification + + start_date = datetime(2018, 1, 1, 0, 0, 0) + end_date = datetime(2018, 12, 31, 0, 0, 0) + timestamp_list = pd.date_range(start_date, end_date, freq='1D') + + # %% Set up logging + logging.basicConfig(level=logging.INFO, + format='%(asctime)s - %(levelname)s - %(message)s', + handlers=[logging.FileHandler(f'ex_opf_{sim_name}.log'), + logging.StreamHandler()]) + + prog_start = time.time() + logging.info(f'Start running multi-period OPF simulation {sim_name}.') + + # %% Set up directories + cwd = os.getcwd() + if 'examples' in cwd: + parent_dir = os.path.dirname(cwd) + data_dir = os.path.join(parent_dir, 'data') + else: + data_dir = os.path.join(cwd, 'data') + + grid_data_dir = os.path.join(data_dir, 'grid') + if not os.path.exists(grid_data_dir): + raise FileNotFoundError('Grid data directory not found.') + + logging.info('Grid data directory: {}'.format(grid_data_dir)) + + fig_dir = os.path.join(os.path.dirname(data_dir), 'figures') + logging.info('Figure directory: {}'.format(fig_dir)) + + results_dir = os.path.join(os.path.dirname(data_dir), 'results') + logging.info('Results directory: {}'.format(results_dir)) + + solar_data_dir = os.path.join(data_dir, 'solar') + print('Solar data directory: {}'.format(solar_data_dir)) + + onshore_wind_data_dir = os.path.join(data_dir, 'onshore_wind') + print('Onshore wind data directory: {}'.format(onshore_wind_data_dir)) + + offshore_wind_data_dir = os.path.join(data_dir, 'offshore_wind') + print('Offshore wind data directory: {}'.format(offshore_wind_data_dir)) + + sim_results_dir = os.path.join(results_dir, sim_name) + if not os.path.exists(sim_results_dir): + os.mkdir(sim_results_dir) + + # %% Read grid data + + # Read load and generation profiles + grid_data = read_grid_data(grid_data_dir, start_date.year) + + # Read DC line property file + filename = os.path.join(grid_data_dir, 'dcline_prop.csv') + dcline_prop = pd.read_csv(filename, index_col=0) + + if w_cpny: + grid_data['dcline_prop'] = dcline_prop + logging.info('With CPNY and CHPE HVDC lines.') # Existing and planned HVDC lines + else: + grid_data['dcline_prop'] = dcline_prop[:4] # Only existing HVDC lines + logging.info('Without CPNY and CHPE HVDC lines.') + + # Read ESR property file + filename = os.path.join(grid_data_dir, 'esr_prop.csv') + esr_prop = pd.read_csv(filename, index_col=0) + + if w_esr: + logging.info('With ESRs.') + grid_data['esr_prop'] = esr_prop # Existing and planned ESRs + else: + logging.info('No ESRs.') + grid_data['esr_prop'] = esr_prop[:8] # Only existing ESRs + + # Read renewable generation profiles + if w_vre: + vre_prop, genmax_profile_vre = read_vre_data(solar_data_dir, onshore_wind_data_dir, offshore_wind_data_dir) + grid_data['vre_prop'] = vre_prop + grid_data['genmax_profile_vre'] = genmax_profile_vre + logging.info('With future solar and offshore wind.') + else: + logging.info('No future solar and offshore wind.') + grid_data['vre_prop'] = None + grid_data['genmax_profile_vre'] = None + + # %% Set up OPF model + + # Set options + options = {'UsePTDF': True, + 'solver': 'gurobi', + 'PenaltyForLoadShed': 20_000, + 'PenaltyForBranchMwViolation': 10_000, + 'PenaltyForInterfaceMWViolation': 10_000} + + # No initial condition for the first day + last_gen = None + + # Loop through all days + for d in range(len(timestamp_list) - 1): + t = time.time() + + # Run OPF for one day (24 hours) plus leading hours + # The first day is valid, the leading hours are used to dispatch batteries properly + start_datetime = timestamp_list[d] + end_datetime = start_datetime + timedelta(hours=24 + leading_hours) + + nygrid_results = run_nygrid_one_day(start_datetime, end_datetime, grid_data, grid_data_dir, options, last_gen) + + # Set generator initial condition for the next iteration + last_gen = nygrid_results['PG'].loc[start_datetime].to_numpy().squeeze() + + # Save simulation nygrid_results to pickle + filename = f'nygrid_sim_{sim_name}_{start_datetime.strftime("%Y%m%d%H")}.pkl' + with open(os.path.join(sim_results_dir, filename), 'wb') as f: + pickle.dump(nygrid_results, f) + logging.info(f'Saved simulation nygrid_results in {filename}') + + elapsed = time.time() - t + logging.info(f'Finished running for {start_datetime.strftime("%Y-%m-%d")}. Elapsed time: {elapsed:.2f} seconds') + logging.info('-' * 80) + + tot_elapsed = time.time() - prog_start + logging.info(f"Finished multi-period OPF simulation {sim_name}. Total elapsed time: {tot_elapsed:.2f} seconds") diff --git a/examples/07_opf_w_vre_cpny.py b/examples/07_opf_w_vre_cpny.py new file mode 100644 index 0000000..67861ae --- /dev/null +++ b/examples/07_opf_w_vre_cpny.py @@ -0,0 +1,153 @@ +""" +Run multi-period OPF with 2018 data. +7. With CPNY and VRE case: +Add CPNY and CHPE HVDC lines. +No ESRs. +Add future solar and offshore wind. +No residential building electrification. +""" +# %% Packages + +import logging +import os +import pickle +import time +from datetime import datetime, timedelta + +import pandas as pd + +from nygrid.run_nygrid import read_grid_data, read_vre_data, run_nygrid_one_day + +if __name__ == '__main__': + + # %% Simulation settings + # NOTE: Change the following settings to run the simulation + sim_name = 'w_vre_cpny' + leading_hours = 12 + w_cpny = True # True: add CPNY and CHPE HVDC lines; False: no CPNY and CHPE HVDC lines + w_esr = False # True: add ESRs; False: no ESRs + w_vre = True # True: add future solar and offshore wind; False: no future solar and offshore wind + w_elec = False # True: add residential building electrification; False: no residential building electrification + + start_date = datetime(2018, 1, 1, 0, 0, 0) + end_date = datetime(2018, 12, 31, 0, 0, 0) + timestamp_list = pd.date_range(start_date, end_date, freq='1D') + + # %% Set up logging + logging.basicConfig(level=logging.INFO, + format='%(asctime)s - %(levelname)s - %(message)s', + handlers=[logging.FileHandler(f'ex_opf_{sim_name}.log'), + logging.StreamHandler()]) + + prog_start = time.time() + logging.info(f'Start running multi-period OPF simulation {sim_name}.') + + # %% Set up directories + cwd = os.getcwd() + if 'examples' in cwd: + parent_dir = os.path.dirname(cwd) + data_dir = os.path.join(parent_dir, 'data') + else: + data_dir = os.path.join(cwd, 'data') + + grid_data_dir = os.path.join(data_dir, 'grid') + if not os.path.exists(grid_data_dir): + raise FileNotFoundError('Grid data directory not found.') + + logging.info('Grid data directory: {}'.format(grid_data_dir)) + + fig_dir = os.path.join(os.path.dirname(data_dir), 'figures') + logging.info('Figure directory: {}'.format(fig_dir)) + + results_dir = os.path.join(os.path.dirname(data_dir), 'results') + logging.info('Results directory: {}'.format(results_dir)) + + solar_data_dir = os.path.join(data_dir, 'solar') + print('Solar data directory: {}'.format(solar_data_dir)) + + onshore_wind_data_dir = os.path.join(data_dir, 'onshore_wind') + print('Onshore wind data directory: {}'.format(onshore_wind_data_dir)) + + offshore_wind_data_dir = os.path.join(data_dir, 'offshore_wind') + print('Offshore wind data directory: {}'.format(offshore_wind_data_dir)) + + sim_results_dir = os.path.join(results_dir, sim_name) + if not os.path.exists(sim_results_dir): + os.mkdir(sim_results_dir) + + # %% Read grid data + + # Read load and generation profiles + grid_data = read_grid_data(grid_data_dir, start_date.year) + + # Read DC line property file + filename = os.path.join(grid_data_dir, 'dcline_prop.csv') + dcline_prop = pd.read_csv(filename, index_col=0) + + if w_cpny: + grid_data['dcline_prop'] = dcline_prop + logging.info('With CPNY and CHPE HVDC lines.') # Existing and planned HVDC lines + else: + grid_data['dcline_prop'] = dcline_prop[:4] # Only existing HVDC lines + logging.info('Without CPNY and CHPE HVDC lines.') + + # Read ESR property file + filename = os.path.join(grid_data_dir, 'esr_prop.csv') + esr_prop = pd.read_csv(filename, index_col=0) + + if w_esr: + logging.info('With ESRs.') + grid_data['esr_prop'] = esr_prop # Existing and planned ESRs + else: + logging.info('No ESRs.') + grid_data['esr_prop'] = esr_prop[:8] # Only existing ESRs + + # Read renewable generation profiles + if w_vre: + vre_prop, genmax_profile_vre = read_vre_data(solar_data_dir, onshore_wind_data_dir, offshore_wind_data_dir) + grid_data['vre_prop'] = vre_prop + grid_data['genmax_profile_vre'] = genmax_profile_vre + logging.info('With future solar and offshore wind.') + else: + logging.info('No future solar and offshore wind.') + grid_data['vre_prop'] = None + grid_data['genmax_profile_vre'] = None + + # %% Set up OPF model + + # Set options + options = {'UsePTDF': True, + 'solver': 'gurobi', + 'PenaltyForLoadShed': 20_000, + 'PenaltyForBranchMwViolation': 10_000, + 'PenaltyForInterfaceMWViolation': 10_000} + + # No initial condition for the first day + last_gen = None + + # Loop through all days + for d in range(len(timestamp_list) - 1): + t = time.time() + + # Run OPF for one day (24 hours) plus leading hours + # The first day is valid, the leading hours are used to dispatch batteries properly + start_datetime = timestamp_list[d] + end_datetime = start_datetime + timedelta(hours=24 + leading_hours) + + nygrid_results = run_nygrid_one_day(start_datetime, end_datetime, grid_data, grid_data_dir, options, last_gen) + + # Set generator initial condition for the next iteration + last_gen = nygrid_results['PG'].loc[start_datetime].to_numpy().squeeze() + + # Save simulation nygrid_results to pickle + filename = f'nygrid_sim_{sim_name}_{start_datetime.strftime("%Y%m%d%H")}.pkl' + with open(os.path.join(sim_results_dir, filename), 'wb') as f: + pickle.dump(nygrid_results, f) + logging.info(f'Saved simulation nygrid_results in {filename}') + + elapsed = time.time() - t + logging.info(f'Finished running for {start_datetime.strftime("%Y-%m-%d")}. Elapsed time: {elapsed:.2f} seconds') + logging.info('-' * 80) + + tot_elapsed = time.time() - prog_start + logging.info(f"Finished multi-period OPF simulation {sim_name}. Total elapsed time: {tot_elapsed:.2f} seconds") diff --git a/examples/08_opf_w_vre_esr_cpny.py b/examples/08_opf_w_vre_esr_cpny.py new file mode 100644 index 0000000..1aa1d37 --- /dev/null +++ b/examples/08_opf_w_vre_esr_cpny.py @@ -0,0 +1,153 @@ +""" +Run multi-period OPF with 2018 data. +8. With CPNY, ESR and VRE case: +Add CPNY and CHPE HVDC lines. +Add 6 GW ESRs. +Add future solar and offshore wind. +No residential building electrification. +""" +# %% Packages + +import logging +import os +import pickle +import time +from datetime import datetime, timedelta + +import pandas as pd + +from nygrid.run_nygrid import read_grid_data, read_vre_data, run_nygrid_one_day + +if __name__ == '__main__': + + # %% Simulation settings + # NOTE: Change the following settings to run the simulation + sim_name = 'w_vre_esr_cpny' + leading_hours = 12 + w_cpny = True # True: add CPNY and CHPE HVDC lines; False: no CPNY and CHPE HVDC lines + w_esr = True # True: add ESRs; False: no ESRs + w_vre = True # True: add future solar and offshore wind; False: no future solar and offshore wind + w_elec = False # True: add residential building electrification; False: no residential building electrification + + start_date = datetime(2018, 1, 1, 0, 0, 0) + end_date = datetime(2018, 12, 31, 0, 0, 0) + timestamp_list = pd.date_range(start_date, end_date, freq='1D') + + # %% Set up logging + logging.basicConfig(level=logging.INFO, + format='%(asctime)s - %(levelname)s - %(message)s', + handlers=[logging.FileHandler(f'ex_opf_{sim_name}.log'), + logging.StreamHandler()]) + + prog_start = time.time() + logging.info(f'Start running multi-period OPF simulation {sim_name}.') + + # %% Set up directories + cwd = os.getcwd() + if 'examples' in cwd: + parent_dir = os.path.dirname(cwd) + data_dir = os.path.join(parent_dir, 'data') + else: + data_dir = os.path.join(cwd, 'data') + + grid_data_dir = os.path.join(data_dir, 'grid') + if not os.path.exists(grid_data_dir): + raise FileNotFoundError('Grid data directory not found.') + + logging.info('Grid data directory: {}'.format(grid_data_dir)) + + fig_dir = os.path.join(os.path.dirname(data_dir), 'figures') + logging.info('Figure directory: {}'.format(fig_dir)) + + results_dir = os.path.join(os.path.dirname(data_dir), 'results') + logging.info('Results directory: {}'.format(results_dir)) + + solar_data_dir = os.path.join(data_dir, 'solar') + print('Solar data directory: {}'.format(solar_data_dir)) + + onshore_wind_data_dir = os.path.join(data_dir, 'onshore_wind') + print('Onshore wind data directory: {}'.format(onshore_wind_data_dir)) + + offshore_wind_data_dir = os.path.join(data_dir, 'offshore_wind') + print('Offshore wind data directory: {}'.format(offshore_wind_data_dir)) + + sim_results_dir = os.path.join(results_dir, sim_name) + if not os.path.exists(sim_results_dir): + os.mkdir(sim_results_dir) + + # %% Read grid data + + # Read load and generation profiles + grid_data = read_grid_data(grid_data_dir, start_date.year) + + # Read DC line property file + filename = os.path.join(grid_data_dir, 'dcline_prop.csv') + dcline_prop = pd.read_csv(filename, index_col=0) + + if w_cpny: + grid_data['dcline_prop'] = dcline_prop + logging.info('With CPNY and CHPE HVDC lines.') # Existing and planned HVDC lines + else: + grid_data['dcline_prop'] = dcline_prop[:4] # Only existing HVDC lines + logging.info('Without CPNY and CHPE HVDC lines.') + + # Read ESR property file + filename = os.path.join(grid_data_dir, 'esr_prop.csv') + esr_prop = pd.read_csv(filename, index_col=0) + + if w_esr: + logging.info('With ESRs.') + grid_data['esr_prop'] = esr_prop # Existing and planned ESRs + else: + logging.info('No ESRs.') + grid_data['esr_prop'] = esr_prop[:8] # Only existing ESRs + + # Read renewable generation profiles + if w_vre: + vre_prop, genmax_profile_vre = read_vre_data(solar_data_dir, onshore_wind_data_dir, offshore_wind_data_dir) + grid_data['vre_prop'] = vre_prop + grid_data['genmax_profile_vre'] = genmax_profile_vre + logging.info('With future solar and offshore wind.') + else: + logging.info('No future solar and offshore wind.') + grid_data['vre_prop'] = None + grid_data['genmax_profile_vre'] = None + + # %% Set up OPF model + + # Set options + options = {'UsePTDF': True, + 'solver': 'gurobi', + 'PenaltyForLoadShed': 20_000, + 'PenaltyForBranchMwViolation': 10_000, + 'PenaltyForInterfaceMWViolation': 10_000} + + # No initial condition for the first day + last_gen = None + + # Loop through all days + for d in range(len(timestamp_list) - 1): + t = time.time() + + # Run OPF for one day (24 hours) plus leading hours + # The first day is valid, the leading hours are used to dispatch batteries properly + start_datetime = timestamp_list[d] + end_datetime = start_datetime + timedelta(hours=24 + leading_hours) + + nygrid_results = run_nygrid_one_day(start_datetime, end_datetime, grid_data, grid_data_dir, options, last_gen) + + # Set generator initial condition for the next iteration + last_gen = nygrid_results['PG'].loc[start_datetime].to_numpy().squeeze() + + # Save simulation nygrid_results to pickle + filename = f'nygrid_sim_{sim_name}_{start_datetime.strftime("%Y%m%d%H")}.pkl' + with open(os.path.join(sim_results_dir, filename), 'wb') as f: + pickle.dump(nygrid_results, f) + logging.info(f'Saved simulation nygrid_results in {filename}') + + elapsed = time.time() - t + logging.info(f'Finished running for {start_datetime.strftime("%Y-%m-%d")}. Elapsed time: {elapsed:.2f} seconds') + logging.info('-' * 80) + + tot_elapsed = time.time() - prog_start + logging.info(f"Finished multi-period OPF simulation {sim_name}. Total elapsed time: {tot_elapsed:.2f} seconds") diff --git a/examples/09_opf_w_elec_vre.py b/examples/09_opf_w_elec_vre.py new file mode 100644 index 0000000..e24d4e1 --- /dev/null +++ b/examples/09_opf_w_elec_vre.py @@ -0,0 +1,169 @@ +""" +Run multi-period OPF with 2018 data. +9. With VRE and electrification case: +No CPNY and CHPE HVDC lines. +No ESRs. +Add future solar and offshore wind. +Add residential building electrification. +""" +# %% Packages + +import logging +import os +import pickle +import time +from datetime import datetime, timedelta + +import pandas as pd + +from nygrid.run_nygrid import read_grid_data, read_vre_data, read_electrification_data, run_nygrid_one_day + +if __name__ == '__main__': + + # %% Simulation settings + # NOTE: Change the following settings to run the simulation + sim_name = 'w_elec_vre' + leading_hours = 12 + w_cpny = False # True: add CPNY and CHPE HVDC lines; False: no CPNY and CHPE HVDC lines + w_esr = False # True: add ESRs; False: no ESRs + w_vre = True # True: add future solar and offshore wind; False: no future solar and offshore wind + w_elec = True # True: add residential building electrification; False: no residential building electrification + + # NOTE: Scale down residential load change by 50% + elec_ratio = 0.5 # Ratio of residential building electrification + + start_date = datetime(2018, 1, 1, 0, 0, 0) + end_date = datetime(2018, 12, 31, 0, 0, 0) + timestamp_list = pd.date_range(start_date, end_date, freq='1D') + + # %% Set up logging + logging.basicConfig(level=logging.INFO, + format='%(asctime)s - %(levelname)s - %(message)s', + handlers=[logging.FileHandler(f'ex_opf_{sim_name}.log'), + logging.StreamHandler()]) + + prog_start = time.time() + logging.info(f'Start running multi-period OPF simulation {sim_name}.') + + # %% Set up directories + cwd = os.getcwd() + if 'examples' in cwd: + parent_dir = os.path.dirname(cwd) + data_dir = os.path.join(parent_dir, 'data') + else: + data_dir = os.path.join(cwd, 'data') + + grid_data_dir = os.path.join(data_dir, 'grid') + if not os.path.exists(grid_data_dir): + raise FileNotFoundError('Grid data directory not found.') + + logging.info('Grid data directory: {}'.format(grid_data_dir)) + + fig_dir = os.path.join(os.path.dirname(data_dir), 'figures') + logging.info('Figure directory: {}'.format(fig_dir)) + + results_dir = os.path.join(os.path.dirname(data_dir), 'results') + logging.info('Results directory: {}'.format(results_dir)) + + solar_data_dir = os.path.join(data_dir, 'solar') + print('Solar data directory: {}'.format(solar_data_dir)) + + onshore_wind_data_dir = os.path.join(data_dir, 'onshore_wind') + print('Onshore wind data directory: {}'.format(onshore_wind_data_dir)) + + offshore_wind_data_dir = os.path.join(data_dir, 'offshore_wind') + print('Offshore wind data directory: {}'.format(offshore_wind_data_dir)) + + buildings_data_dir = os.path.join(data_dir, 'buildings') + print('Buildings data directory: {}'.format(buildings_data_dir)) + + sim_results_dir = os.path.join(results_dir, sim_name) + if not os.path.exists(sim_results_dir): + os.mkdir(sim_results_dir) + + # %% Read grid data + + # Read load and generation profiles + grid_data = read_grid_data(grid_data_dir, start_date.year) + + # Read DC line property file + filename = os.path.join(grid_data_dir, 'dcline_prop.csv') + dcline_prop = pd.read_csv(filename, index_col=0) + + if w_cpny: + grid_data['dcline_prop'] = dcline_prop + logging.info('With CPNY and CHPE HVDC lines.') # Existing and planned HVDC lines + else: + grid_data['dcline_prop'] = dcline_prop[:4] # Only existing HVDC lines + logging.info('Without CPNY and CHPE HVDC lines.') + + # Read ESR property file + filename = os.path.join(grid_data_dir, 'esr_prop.csv') + esr_prop = pd.read_csv(filename, index_col=0) + + if w_esr: + logging.info('With ESRs.') + grid_data['esr_prop'] = esr_prop # Existing and planned ESRs + else: + logging.info('No ESRs.') + grid_data['esr_prop'] = esr_prop[:8] # Only existing ESRs + + # Read renewable generation profiles + if w_vre: + vre_prop, genmax_profile_vre = read_vre_data(solar_data_dir, onshore_wind_data_dir, offshore_wind_data_dir) + grid_data['vre_prop'] = vre_prop + grid_data['genmax_profile_vre'] = genmax_profile_vre + logging.info('With future solar and offshore wind.') + else: + logging.info('No future solar and offshore wind.') + grid_data['vre_prop'] = None + grid_data['genmax_profile_vre'] = None + + # Read residential building electrification profiles + if w_elec: + logging.info('With residential building electrification.') + res_load_change_bus = read_electrification_data(buildings_data_dir) + res_load_change_bus = res_load_change_bus * elec_ratio + load_profile_elec = grid_data['load_profile'].add(res_load_change_bus, fill_value=0) + grid_data['load_profile'] = load_profile_elec + else: + logging.info('No residential building electrification.') + + # %% Set up OPF model + + # Set options + options = {'UsePTDF': True, + 'solver': 'gurobi', + 'PenaltyForLoadShed': 20_000, + 'PenaltyForBranchMwViolation': 10_000, + 'PenaltyForInterfaceMWViolation': 10_000} + + # No initial condition for the first day + last_gen = None + + # Loop through all days + for d in range(len(timestamp_list) - 1): + t = time.time() + + # Run OPF for one day (24 hours) plus leading hours + # The first day is valid, the leading hours are used to dispatch batteries properly + start_datetime = timestamp_list[d] + end_datetime = start_datetime + timedelta(hours=24 + leading_hours) + + nygrid_results = run_nygrid_one_day(start_datetime, end_datetime, grid_data, grid_data_dir, options, last_gen) + + # Set generator initial condition for the next iteration + last_gen = nygrid_results['PG'].loc[start_datetime].to_numpy().squeeze() + + # Save simulation nygrid_results to pickle + filename = f'nygrid_sim_{sim_name}_{start_datetime.strftime("%Y%m%d%H")}.pkl' + with open(os.path.join(sim_results_dir, filename), 'wb') as f: + pickle.dump(nygrid_results, f) + logging.info(f'Saved simulation nygrid_results in {filename}') + + elapsed = time.time() - t + logging.info(f'Finished running for {start_datetime.strftime("%Y-%m-%d")}. Elapsed time: {elapsed:.2f} seconds') + logging.info('-' * 80) + + tot_elapsed = time.time() - prog_start + logging.info(f"Finished multi-period OPF simulation {sim_name}. Total elapsed time: {tot_elapsed:.2f} seconds") diff --git a/examples/10_opf_w_elec_vre_esr.py b/examples/10_opf_w_elec_vre_esr.py new file mode 100644 index 0000000..d621d17 --- /dev/null +++ b/examples/10_opf_w_elec_vre_esr.py @@ -0,0 +1,169 @@ +""" +Run multi-period OPF with 2018 data. +10. With ESR, VRE and electrification case: +No CPNY and CHPE HVDC lines. +Add 6 GW ESRs. +Add future solar and offshore wind. +Add residential building electrification. +""" +# %% Packages + +import logging +import os +import pickle +import time +from datetime import datetime, timedelta + +import pandas as pd + +from nygrid.run_nygrid import read_grid_data, read_vre_data, read_electrification_data, run_nygrid_one_day + +if __name__ == '__main__': + + # %% Simulation settings + # NOTE: Change the following settings to run the simulation + sim_name = 'w_elec_vre_esr' + leading_hours = 12 + w_cpny = False # True: add CPNY and CHPE HVDC lines; False: no CPNY and CHPE HVDC lines + w_esr = True # True: add ESRs; False: no ESRs + w_vre = True # True: add future solar and offshore wind; False: no future solar and offshore wind + w_elec = True # True: add residential building electrification; False: no residential building electrification + + # NOTE: Scale down residential load change by 50% + elec_ratio = 0.5 # Ratio of residential building electrification + + start_date = datetime(2018, 1, 1, 0, 0, 0) + end_date = datetime(2018, 12, 31, 0, 0, 0) + timestamp_list = pd.date_range(start_date, end_date, freq='1D') + + # %% Set up logging + logging.basicConfig(level=logging.INFO, + format='%(asctime)s - %(levelname)s - %(message)s', + handlers=[logging.FileHandler(f'ex_opf_{sim_name}.log'), + logging.StreamHandler()]) + + prog_start = time.time() + logging.info(f'Start running multi-period OPF simulation {sim_name}.') + + # %% Set up directories + cwd = os.getcwd() + if 'examples' in cwd: + parent_dir = os.path.dirname(cwd) + data_dir = os.path.join(parent_dir, 'data') + else: + data_dir = os.path.join(cwd, 'data') + + grid_data_dir = os.path.join(data_dir, 'grid') + if not os.path.exists(grid_data_dir): + raise FileNotFoundError('Grid data directory not found.') + + logging.info('Grid data directory: {}'.format(grid_data_dir)) + + fig_dir = os.path.join(os.path.dirname(data_dir), 'figures') + logging.info('Figure directory: {}'.format(fig_dir)) + + results_dir = os.path.join(os.path.dirname(data_dir), 'results') + logging.info('Results directory: {}'.format(results_dir)) + + solar_data_dir = os.path.join(data_dir, 'solar') + print('Solar data directory: {}'.format(solar_data_dir)) + + onshore_wind_data_dir = os.path.join(data_dir, 'onshore_wind') + print('Onshore wind data directory: {}'.format(onshore_wind_data_dir)) + + offshore_wind_data_dir = os.path.join(data_dir, 'offshore_wind') + print('Offshore wind data directory: {}'.format(offshore_wind_data_dir)) + + buildings_data_dir = os.path.join(data_dir, 'buildings') + print('Buildings data directory: {}'.format(buildings_data_dir)) + + sim_results_dir = os.path.join(results_dir, sim_name) + if not os.path.exists(sim_results_dir): + os.mkdir(sim_results_dir) + + # %% Read grid data + + # Read load and generation profiles + grid_data = read_grid_data(grid_data_dir, start_date.year) + + # Read DC line property file + filename = os.path.join(grid_data_dir, 'dcline_prop.csv') + dcline_prop = pd.read_csv(filename, index_col=0) + + if w_cpny: + grid_data['dcline_prop'] = dcline_prop + logging.info('With CPNY and CHPE HVDC lines.') # Existing and planned HVDC lines + else: + grid_data['dcline_prop'] = dcline_prop[:4] # Only existing HVDC lines + logging.info('Without CPNY and CHPE HVDC lines.') + + # Read ESR property file + filename = os.path.join(grid_data_dir, 'esr_prop.csv') + esr_prop = pd.read_csv(filename, index_col=0) + + if w_esr: + logging.info('With ESRs.') + grid_data['esr_prop'] = esr_prop # Existing and planned ESRs + else: + logging.info('No ESRs.') + grid_data['esr_prop'] = esr_prop[:8] # Only existing ESRs + + # Read renewable generation profiles + if w_vre: + vre_prop, genmax_profile_vre = read_vre_data(solar_data_dir, onshore_wind_data_dir, offshore_wind_data_dir) + grid_data['vre_prop'] = vre_prop + grid_data['genmax_profile_vre'] = genmax_profile_vre + logging.info('With future solar and offshore wind.') + else: + logging.info('No future solar and offshore wind.') + grid_data['vre_prop'] = None + grid_data['genmax_profile_vre'] = None + + # Read residential building electrification profiles + if w_elec: + logging.info('With residential building electrification.') + res_load_change_bus = read_electrification_data(buildings_data_dir) + res_load_change_bus = res_load_change_bus * elec_ratio + load_profile_elec = grid_data['load_profile'].add(res_load_change_bus, fill_value=0) + grid_data['load_profile'] = load_profile_elec + else: + logging.info('No residential building electrification.') + + # %% Set up OPF model + + # Set options + options = {'UsePTDF': True, + 'solver': 'gurobi', + 'PenaltyForLoadShed': 20_000, + 'PenaltyForBranchMwViolation': 10_000, + 'PenaltyForInterfaceMWViolation': 10_000} + + # No initial condition for the first day + last_gen = None + + # Loop through all days + for d in range(len(timestamp_list) - 1): + t = time.time() + + # Run OPF for one day (24 hours) plus leading hours + # The first day is valid, the leading hours are used to dispatch batteries properly + start_datetime = timestamp_list[d] + end_datetime = start_datetime + timedelta(hours=24 + leading_hours) + + nygrid_results = run_nygrid_one_day(start_datetime, end_datetime, grid_data, grid_data_dir, options, last_gen) + + # Set generator initial condition for the next iteration + last_gen = nygrid_results['PG'].loc[start_datetime].to_numpy().squeeze() + + # Save simulation nygrid_results to pickle + filename = f'nygrid_sim_{sim_name}_{start_datetime.strftime("%Y%m%d%H")}.pkl' + with open(os.path.join(sim_results_dir, filename), 'wb') as f: + pickle.dump(nygrid_results, f) + logging.info(f'Saved simulation nygrid_results in {filename}') + + elapsed = time.time() - t + logging.info(f'Finished running for {start_datetime.strftime("%Y-%m-%d")}. Elapsed time: {elapsed:.2f} seconds') + logging.info('-' * 80) + + tot_elapsed = time.time() - prog_start + logging.info(f"Finished multi-period OPF simulation {sim_name}. Total elapsed time: {tot_elapsed:.2f} seconds") diff --git a/examples/11_opf_w_elec_vre_cpny.py b/examples/11_opf_w_elec_vre_cpny.py new file mode 100644 index 0000000..725a231 --- /dev/null +++ b/examples/11_opf_w_elec_vre_cpny.py @@ -0,0 +1,169 @@ +""" +Run multi-period OPF with 2018 data. +11. With CPNY, VRE and electrification case: +Add CPNY and CHPE HVDC lines. +No ESRs. +Add future solar and offshore wind. +Add residential building electrification. +""" +# %% Packages + +import logging +import os +import pickle +import time +from datetime import datetime, timedelta + +import pandas as pd + +from nygrid.run_nygrid import read_grid_data, read_vre_data, read_electrification_data, run_nygrid_one_day + +if __name__ == '__main__': + + # %% Simulation settings + # NOTE: Change the following settings to run the simulation + sim_name = 'w_elec_vre_cpny' + leading_hours = 12 + w_cpny = True # True: add CPNY and CHPE HVDC lines; False: no CPNY and CHPE HVDC lines + w_esr = False # True: add ESRs; False: no ESRs + w_vre = True # True: add future solar and offshore wind; False: no future solar and offshore wind + w_elec = True # True: add residential building electrification; False: no residential building electrification + + # NOTE: Scale down residential load change by 50% + elec_ratio = 0.5 # Ratio of residential building electrification + + start_date = datetime(2018, 1, 1, 0, 0, 0) + end_date = datetime(2018, 12, 31, 0, 0, 0) + timestamp_list = pd.date_range(start_date, end_date, freq='1D') + + # %% Set up logging + logging.basicConfig(level=logging.INFO, + format='%(asctime)s - %(levelname)s - %(message)s', + handlers=[logging.FileHandler(f'ex_opf_{sim_name}.log'), + logging.StreamHandler()]) + + prog_start = time.time() + logging.info(f'Start running multi-period OPF simulation {sim_name}.') + + # %% Set up directories + cwd = os.getcwd() + if 'examples' in cwd: + parent_dir = os.path.dirname(cwd) + data_dir = os.path.join(parent_dir, 'data') + else: + data_dir = os.path.join(cwd, 'data') + + grid_data_dir = os.path.join(data_dir, 'grid') + if not os.path.exists(grid_data_dir): + raise FileNotFoundError('Grid data directory not found.') + + logging.info('Grid data directory: {}'.format(grid_data_dir)) + + fig_dir = os.path.join(os.path.dirname(data_dir), 'figures') + logging.info('Figure directory: {}'.format(fig_dir)) + + results_dir = os.path.join(os.path.dirname(data_dir), 'results') + logging.info('Results directory: {}'.format(results_dir)) + + solar_data_dir = os.path.join(data_dir, 'solar') + print('Solar data directory: {}'.format(solar_data_dir)) + + onshore_wind_data_dir = os.path.join(data_dir, 'onshore_wind') + print('Onshore wind data directory: {}'.format(onshore_wind_data_dir)) + + offshore_wind_data_dir = os.path.join(data_dir, 'offshore_wind') + print('Offshore wind data directory: {}'.format(offshore_wind_data_dir)) + + buildings_data_dir = os.path.join(data_dir, 'buildings') + print('Buildings data directory: {}'.format(buildings_data_dir)) + + sim_results_dir = os.path.join(results_dir, sim_name) + if not os.path.exists(sim_results_dir): + os.mkdir(sim_results_dir) + + # %% Read grid data + + # Read load and generation profiles + grid_data = read_grid_data(grid_data_dir, start_date.year) + + # Read DC line property file + filename = os.path.join(grid_data_dir, 'dcline_prop.csv') + dcline_prop = pd.read_csv(filename, index_col=0) + + if w_cpny: + grid_data['dcline_prop'] = dcline_prop + logging.info('With CPNY and CHPE HVDC lines.') # Existing and planned HVDC lines + else: + grid_data['dcline_prop'] = dcline_prop[:4] # Only existing HVDC lines + logging.info('Without CPNY and CHPE HVDC lines.') + + # Read ESR property file + filename = os.path.join(grid_data_dir, 'esr_prop.csv') + esr_prop = pd.read_csv(filename, index_col=0) + + if w_esr: + logging.info('With ESRs.') + grid_data['esr_prop'] = esr_prop # Existing and planned ESRs + else: + logging.info('No ESRs.') + grid_data['esr_prop'] = esr_prop[:8] # Only existing ESRs + + # Read renewable generation profiles + if w_vre: + vre_prop, genmax_profile_vre = read_vre_data(solar_data_dir, onshore_wind_data_dir, offshore_wind_data_dir) + grid_data['vre_prop'] = vre_prop + grid_data['genmax_profile_vre'] = genmax_profile_vre + logging.info('With future solar and offshore wind.') + else: + logging.info('No future solar and offshore wind.') + grid_data['vre_prop'] = None + grid_data['genmax_profile_vre'] = None + + # Read residential building electrification profiles + if w_elec: + logging.info('With residential building electrification.') + res_load_change_bus = read_electrification_data(buildings_data_dir) + res_load_change_bus = res_load_change_bus * elec_ratio + load_profile_elec = grid_data['load_profile'].add(res_load_change_bus, fill_value=0) + grid_data['load_profile'] = load_profile_elec + else: + logging.info('No residential building electrification.') + + # %% Set up OPF model + + # Set options + options = {'UsePTDF': True, + 'solver': 'gurobi', + 'PenaltyForLoadShed': 20_000, + 'PenaltyForBranchMwViolation': 10_000, + 'PenaltyForInterfaceMWViolation': 10_000} + + # No initial condition for the first day + last_gen = None + + # Loop through all days + for d in range(len(timestamp_list) - 1): + t = time.time() + + # Run OPF for one day (24 hours) plus leading hours + # The first day is valid, the leading hours are used to dispatch batteries properly + start_datetime = timestamp_list[d] + end_datetime = start_datetime + timedelta(hours=24 + leading_hours) + + nygrid_results = run_nygrid_one_day(start_datetime, end_datetime, grid_data, grid_data_dir, options, last_gen) + + # Set generator initial condition for the next iteration + last_gen = nygrid_results['PG'].loc[start_datetime].to_numpy().squeeze() + + # Save simulation nygrid_results to pickle + filename = f'nygrid_sim_{sim_name}_{start_datetime.strftime("%Y%m%d%H")}.pkl' + with open(os.path.join(sim_results_dir, filename), 'wb') as f: + pickle.dump(nygrid_results, f) + logging.info(f'Saved simulation nygrid_results in {filename}') + + elapsed = time.time() - t + logging.info(f'Finished running for {start_datetime.strftime("%Y-%m-%d")}. Elapsed time: {elapsed:.2f} seconds') + logging.info('-' * 80) + + tot_elapsed = time.time() - prog_start + logging.info(f"Finished multi-period OPF simulation {sim_name}. Total elapsed time: {tot_elapsed:.2f} seconds") diff --git a/examples/12_opf_w_elec_vre_esr_cpny.py b/examples/12_opf_w_elec_vre_esr_cpny.py new file mode 100644 index 0000000..92c3e6e --- /dev/null +++ b/examples/12_opf_w_elec_vre_esr_cpny.py @@ -0,0 +1,169 @@ +""" +Run multi-period OPF with 2018 data. +12. With ESR, CPNY, VRE and electrification case: +Add CPNY and CHPE HVDC lines. +Add 6 GW ESRs. +Add future solar and offshore wind. +Add residential building electrification. +""" +# %% Packages + +import logging +import os +import pickle +import time +from datetime import datetime, timedelta + +import pandas as pd + +from nygrid.run_nygrid import read_grid_data, read_vre_data, read_electrification_data, run_nygrid_one_day + +if __name__ == '__main__': + + # %% Simulation settings + # NOTE: Change the following settings to run the simulation + sim_name = 'w_elec_vre_esr_cpny' + leading_hours = 12 + w_cpny = True # True: add CPNY and CHPE HVDC lines; False: no CPNY and CHPE HVDC lines + w_esr = True # True: add ESRs; False: no ESRs + w_vre = True # True: add future solar and offshore wind; False: no future solar and offshore wind + w_elec = True # True: add residential building electrification; False: no residential building electrification + + # NOTE: Scale down residential load change by 50% + elec_ratio = 0.5 # Ratio of residential building electrification + + start_date = datetime(2018, 1, 1, 0, 0, 0) + end_date = datetime(2018, 12, 31, 0, 0, 0) + timestamp_list = pd.date_range(start_date, end_date, freq='1D') + + # %% Set up logging + logging.basicConfig(level=logging.INFO, + format='%(asctime)s - %(levelname)s - %(message)s', + handlers=[logging.FileHandler(f'ex_opf_{sim_name}.log'), + logging.StreamHandler()]) + + prog_start = time.time() + logging.info(f'Start running multi-period OPF simulation {sim_name}.') + + # %% Set up directories + cwd = os.getcwd() + if 'examples' in cwd: + parent_dir = os.path.dirname(cwd) + data_dir = os.path.join(parent_dir, 'data') + else: + data_dir = os.path.join(cwd, 'data') + + grid_data_dir = os.path.join(data_dir, 'grid') + if not os.path.exists(grid_data_dir): + raise FileNotFoundError('Grid data directory not found.') + + logging.info('Grid data directory: {}'.format(grid_data_dir)) + + fig_dir = os.path.join(os.path.dirname(data_dir), 'figures') + logging.info('Figure directory: {}'.format(fig_dir)) + + results_dir = os.path.join(os.path.dirname(data_dir), 'results') + logging.info('Results directory: {}'.format(results_dir)) + + solar_data_dir = os.path.join(data_dir, 'solar') + print('Solar data directory: {}'.format(solar_data_dir)) + + onshore_wind_data_dir = os.path.join(data_dir, 'onshore_wind') + print('Onshore wind data directory: {}'.format(onshore_wind_data_dir)) + + offshore_wind_data_dir = os.path.join(data_dir, 'offshore_wind') + print('Offshore wind data directory: {}'.format(offshore_wind_data_dir)) + + buildings_data_dir = os.path.join(data_dir, 'buildings') + print('Buildings data directory: {}'.format(buildings_data_dir)) + + sim_results_dir = os.path.join(results_dir, sim_name) + if not os.path.exists(sim_results_dir): + os.mkdir(sim_results_dir) + + # %% Read grid data + + # Read load and generation profiles + grid_data = read_grid_data(grid_data_dir, start_date.year) + + # Read DC line property file + filename = os.path.join(grid_data_dir, 'dcline_prop.csv') + dcline_prop = pd.read_csv(filename, index_col=0) + + if w_cpny: + grid_data['dcline_prop'] = dcline_prop + logging.info('With CPNY and CHPE HVDC lines.') # Existing and planned HVDC lines + else: + grid_data['dcline_prop'] = dcline_prop[:4] # Only existing HVDC lines + logging.info('Without CPNY and CHPE HVDC lines.') + + # Read ESR property file + filename = os.path.join(grid_data_dir, 'esr_prop.csv') + esr_prop = pd.read_csv(filename, index_col=0) + + if w_esr: + logging.info('With ESRs.') + grid_data['esr_prop'] = esr_prop # Existing and planned ESRs + else: + logging.info('No ESRs.') + grid_data['esr_prop'] = esr_prop[:8] # Only existing ESRs + + # Read renewable generation profiles + if w_vre: + vre_prop, genmax_profile_vre = read_vre_data(solar_data_dir, onshore_wind_data_dir, offshore_wind_data_dir) + grid_data['vre_prop'] = vre_prop + grid_data['genmax_profile_vre'] = genmax_profile_vre + logging.info('With future solar and offshore wind.') + else: + logging.info('No future solar and offshore wind.') + grid_data['vre_prop'] = None + grid_data['genmax_profile_vre'] = None + + # Read residential building electrification profiles + if w_elec: + logging.info('With residential building electrification.') + res_load_change_bus = read_electrification_data(buildings_data_dir) + res_load_change_bus = res_load_change_bus * elec_ratio + load_profile_elec = grid_data['load_profile'].add(res_load_change_bus, fill_value=0) + grid_data['load_profile'] = load_profile_elec + else: + logging.info('No residential building electrification.') + + # %% Set up OPF model + + # Set options + options = {'UsePTDF': True, + 'solver': 'gurobi', + 'PenaltyForLoadShed': 20_000, + 'PenaltyForBranchMwViolation': 10_000, + 'PenaltyForInterfaceMWViolation': 10_000} + + # No initial condition for the first day + last_gen = None + + # Loop through all days + for d in range(len(timestamp_list) - 1): + t = time.time() + + # Run OPF for one day (24 hours) plus leading hours + # The first day is valid, the leading hours are used to dispatch batteries properly + start_datetime = timestamp_list[d] + end_datetime = start_datetime + timedelta(hours=24 + leading_hours) + + nygrid_results = run_nygrid_one_day(start_datetime, end_datetime, grid_data, grid_data_dir, options, last_gen) + + # Set generator initial condition for the next iteration + last_gen = nygrid_results['PG'].loc[start_datetime].to_numpy().squeeze() + + # Save simulation nygrid_results to pickle + filename = f'nygrid_sim_{sim_name}_{start_datetime.strftime("%Y%m%d%H")}.pkl' + with open(os.path.join(sim_results_dir, filename), 'wb') as f: + pickle.dump(nygrid_results, f) + logging.info(f'Saved simulation nygrid_results in {filename}') + + elapsed = time.time() - t + logging.info(f'Finished running for {start_datetime.strftime("%Y-%m-%d")}. Elapsed time: {elapsed:.2f} seconds') + logging.info('-' * 80) + + tot_elapsed = time.time() - prog_start + logging.info(f"Finished multi-period OPF simulation {sim_name}. Total elapsed time: {tot_elapsed:.2f} seconds") diff --git a/examples/ex_opf_w_future_solar.py b/examples/ex_opf_w_future_solar.py deleted file mode 100644 index b590148..0000000 --- a/examples/ex_opf_w_future_solar.py +++ /dev/null @@ -1,290 +0,0 @@ -""" -Run multi-period OPF with 2018 data -with renewable generators -including wind and solar - -""" -# %% Packages -import os -import sys -import numpy as np -import pandas as pd -from datetime import datetime, timedelta - -from pyomo.opt import SolverFactory - -from nygrid.nygrid import NYGrid, check_status -import pickle -import time -import logging - - -# %% Functions - -def read_renewable_data(solar_data_dir, onshore_wind_data_dir, offshore_wind_data_dir): - # Renewable generation time series - current_solar_gen = pd.read_csv(os.path.join(solar_data_dir, f'current_solar_gen_1hr.csv'), - parse_dates=['Time'], index_col='Time') - current_solar_gen.index.freq = 'H' - current_solar_gen.columns = current_solar_gen.columns.astype(int) - future_solar_gen = pd.read_csv(os.path.join(solar_data_dir, f'future_solar_gen_1hr.csv'), - parse_dates=['Time'], index_col='Time') - future_solar_gen.index.freq = 'H' - future_solar_gen.columns = future_solar_gen.columns.astype(int) - onshore_wind_gen = pd.read_csv(os.path.join(onshore_wind_data_dir, f'current_wind_gen_1hr.csv'), - parse_dates=['Time'], index_col='Time') - onshore_wind_gen.index.freq = 'H' - onshore_wind_gen.columns = onshore_wind_gen.columns.astype(int) - offshore_wind_gen = pd.read_csv(os.path.join(offshore_wind_data_dir, f'power_load_2018.csv'), - parse_dates=['timestamp'], index_col='timestamp') - offshore_wind_gen.index = offshore_wind_gen.index.tz_localize('US/Eastern', ambiguous='infer') - offshore_wind_gen.index.freq = 'H' - # Renewable allocation table - current_solar_2bus = pd.read_csv(os.path.join(solar_data_dir, f'solar_farms_2bus.csv'), - index_col='zip_code') - future_solar_2bus = pd.read_csv(os.path.join(solar_data_dir, f'future_solar_farms_2bus.csv'), index_col=0) - onshore_wind_2bus = pd.read_csv(os.path.join(onshore_wind_data_dir, f'onshore_wind_2bus.csv'), index_col=0) - # Aggregate current solar generation - groupby_dict = current_solar_2bus['busIdx'].to_dict() - current_solar_gen_agg = current_solar_gen.groupby(groupby_dict, axis=1).sum() - current_solar_gen_agg = current_solar_gen_agg / 1e3 # convert from kW to MW - # Aggregate future solar generation - groupby_dict = future_solar_2bus['busIdx'].to_dict() - future_solar_gen_agg = future_solar_gen.groupby(groupby_dict, axis=1).sum() - future_solar_gen_agg = future_solar_gen_agg / 1e3 # convert from kW to MW - # Aggregate onshore wind generation - groupby_dict = onshore_wind_2bus['busIdx'].to_dict() - onshore_wind_gen_agg = onshore_wind_gen.groupby(groupby_dict, axis=1).sum() - onshore_wind_gen_agg = onshore_wind_gen_agg / 1e3 # convert from kW to MW - # Aggregate offshore wind generation - offshore_wind_gen_agg = offshore_wind_gen[['power_nyc', 'power_li']].rename( - columns={'power_nyc': 81, 'power_li': 79}) - - return current_solar_gen_agg, onshore_wind_gen_agg, offshore_wind_gen_agg, future_solar_gen_agg - - -def update_load_profile(load_profile, current_solar_gen_agg=None, onshore_wind_gen_agg=None, future_solar_gen_agg=None): - # Treat renewable generation as negative load - load_profile_copy = load_profile.copy() - load_profile_copy.columns = [int(col.replace('Bus', '')) for col in load_profile_copy.columns] - load_profile_copy.index = current_solar_gen_agg.index - # 18.08% of current solar generation is built before 2018 (base year) - # Scale down current solar generation by 18.08% - pct_current_solar_built = 0.1808 - current_solar_gen_agg = current_solar_gen_agg * (1 - pct_current_solar_built) - # 90.6% of onshore wind generation is built before 2018 (base year) - # Scale down onshore wind generation by 90.6% - pct_onshore_wind_built = 0.906 - onshore_wind_gen_agg = onshore_wind_gen_agg * (1 - pct_onshore_wind_built) - # Total renewable generation by bus - total_renewable = pd.DataFrame() - # total_renewable = total_renewable.add(onshore_wind_gen_agg, fill_value=0) - # total_renewable = total_renewable.add(offshore_wind_gen_agg, fill_value=0) - total_renewable = total_renewable.add(current_solar_gen_agg, fill_value=0) - total_renewable = total_renewable.add(future_solar_gen_agg, fill_value=0) - # Load profile after subtracting renewable generation - load_profile_renewable = load_profile_copy.subtract(total_renewable, fill_value=0) - # Scale down external load by the same percentage as NYS - ny_buses = np.arange(37, 83) - ext_buses = np.array(list(set(load_profile_copy.columns) - set(ny_buses))) - ny_load_profile = load_profile_copy[ny_buses] - ext_load_profile = load_profile_copy[ext_buses] - # NY load change rate in each hour - ny_change_pct = ((load_profile_renewable[ny_buses].sum(axis=1) - - ny_load_profile.sum(axis=1)) / ny_load_profile.sum(axis=1)) - print(f'NYS annual load change: {ny_change_pct.mean() * 100:.2f}%') - ext_change = ext_load_profile.multiply(ny_change_pct * 0.2, axis=0) - load_profile_renewable = load_profile_renewable.add(ext_change, fill_value=0) - print(f'External load change is set to 20% of NYS load change.') - # Reset timestamp index to time zone unaware - load_profile_copy.index = load_profile.index - load_profile_renewable.index = load_profile.index - - return load_profile_renewable - - -def read_grid_data(grid_data_dir, start_date): - # Read load profile - load_profile = pd.read_csv(os.path.join(grid_data_dir, f'load_profile_{start_date.year}.csv'), - parse_dates=['TimeStamp'], index_col='TimeStamp') - load_profile.index.freq = 'H' - - # Read generation profile - gen_profile = pd.read_csv(os.path.join(grid_data_dir, f'gen_profile_{start_date.year}.csv'), - parse_dates=['TimeStamp'], index_col='TimeStamp') - gen_profile.index.freq = 'H' - - # Read generator capacity limit profile - gen_max_profile = pd.read_csv(os.path.join(grid_data_dir, f'genmax_profile_{start_date.year}.csv'), - parse_dates=['TimeStamp'], index_col='TimeStamp') - gen_max_profile.index.freq = 'H' - - gen_min_profile = pd.read_csv(os.path.join(grid_data_dir, f'genmin_profile_{start_date.year}.csv'), - parse_dates=['TimeStamp'], index_col='TimeStamp') - gen_min_profile.index.freq = 'H' - - # Read generator ramp rate profile - gen_ramp30_profile = pd.read_csv(os.path.join(grid_data_dir, f'genramp30_profile_{start_date.year}.csv'), - parse_dates=['TimeStamp'], index_col='TimeStamp') - gen_ramp30_profile.index.freq = 'H' - - # Read generator cost profile (linear) - gen_cost0_profile = pd.read_csv(os.path.join(grid_data_dir, f'gencost0_profile_{start_date.year}.csv'), - parse_dates=['TimeStamp'], index_col='TimeStamp') - gen_cost0_profile.index.freq = 'H' - - gen_cost1_profile = pd.read_csv(os.path.join(grid_data_dir, f'gencost1_profile_{start_date.year}.csv'), - parse_dates=['TimeStamp'], index_col='TimeStamp') - gen_cost1_profile.index.freq = 'H' - - return (load_profile, gen_profile, gen_max_profile, gen_min_profile, - gen_ramp30_profile, gen_cost0_profile, gen_cost1_profile) - - -if __name__ == '__main__': - - # %% Set up the program - - # Set up logging - logging.basicConfig(stream=sys.stdout, level=logging.INFO, - format='%(asctime)s - %(levelname)s - %(message)s') - - # Start timer - t = time.time() - logging.info(f"Starting program at {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}") - - # Set up directories - cwd = os.getcwd() - if 'examples' in cwd: - parent_dir = os.path.dirname(cwd) - data_dir = os.path.join(parent_dir, 'data') - else: - data_dir = os.path.join(cwd, 'data') - - grid_data_dir = os.path.join(data_dir, 'grid') - if not os.path.exists(grid_data_dir): - raise FileNotFoundError('Grid data directory not found.') - - logging.info('Grid data directory: {}'.format(grid_data_dir)) - - fig_dir = os.path.join(os.path.dirname(data_dir), 'figures') - logging.info('Figure directory: {}'.format(fig_dir)) - - results_dir = os.path.join(os.path.dirname(data_dir), 'results') - logging.info('Results directory: {}'.format(results_dir)) - - if not os.path.exists(os.path.join(results_dir, 'w_solar')): - os.mkdir(os.path.join(results_dir, 'w_solar')) - - solar_data_dir = os.path.join(data_dir, 'solar') - logging.info('Solar data directory: {}'.format(solar_data_dir)) - - onshore_wind_data_dir = os.path.join(data_dir, 'onshore_wind') - logging.info('Onshore wind data directory: {}'.format(onshore_wind_data_dir)) - - offshore_wind_data_dir = os.path.join(data_dir, 'offshore_wind') - logging.info('Offshore wind data directory: {}'.format(offshore_wind_data_dir)) - - # %% Read grid data - start_date = datetime(2018, 1, 1, 0, 0, 0) - end_date = datetime(2019, 1, 1, 0, 0, 0) - timestamp_list = pd.date_range(start_date, end_date, freq='1D') - - load_profile, gen_profile, genmax_profile, genmin_profile, genramp30_profile, \ - gencost0_profile, gencost1_profile = read_grid_data(grid_data_dir, start_date) - - # %% Read renewable generation data - - current_solar_gen_agg, onshore_wind_gen_agg, offshore_wind_gen_agg, future_solar_gen_agg = \ - read_renewable_data(solar_data_dir, onshore_wind_data_dir, offshore_wind_data_dir) - - # %% Update load profile - - load_profile_renewable = update_load_profile(load_profile, current_solar_gen_agg, - onshore_wind_gen_agg, future_solar_gen_agg) - - # %% Read thermal generator info table - - filename = os.path.join(data_dir, 'genInfo.csv') - gen_info = pd.read_csv(filename) - num_thermal = gen_info.shape[0] - gen_rename = {gen_info.index[i]: gen_info.NYISOName[i] for i in range(num_thermal)} - - # %% Set up OPF model - timestamp_list = pd.date_range(start_date, end_date, freq='1D') - - # Loop through all days - for d in range(len(timestamp_list) - 1): - - # Run one day at a time - start_datetime = timestamp_list[d] - end_datetime = start_datetime + timedelta(hours=23) - - # Read MATPOWER case file - ppc_filename = os.path.join(data_dir, 'ny_grid.mat') - - nygrid_sim = NYGrid(ppc_filename, - start_datetime=start_datetime.strftime('%m-%d-%Y %H'), - end_datetime=end_datetime.strftime('%m-%d-%Y %H'), - verbose=True) - - # Read grid data - nygrid_sim.set_load_sch(load_profile_renewable) - nygrid_sim.set_gen_mw_sch(gen_profile) - nygrid_sim.set_gen_max_sch(genmax_profile) - nygrid_sim.set_gen_min_sch(genmin_profile) - nygrid_sim.set_gen_ramp_sch(genramp30_profile) - nygrid_sim.set_gen_cost_sch(gencost0_profile, gencost1_profile) - - # Process ppc - nygrid_sim.process_ppc() - - # Set generator initial condition - if d == 0: - last_gen = None - logging.info('No initial condition.') - else: - last_gen = nygrid_sim.get_last_gen(model_multi_opf) - logging.info('Initial condition set from previous day.') - - nygrid_sim.set_gen_init_data(gen_init=last_gen) - - # Check input - nygrid_sim.check_input_dim() - - # Initialize single period OPF - model_multi_opf = nygrid_sim.create_multi_opf_soft(slack_cost_weight=1e21) - - solver = SolverFactory('gurobi') - results_multi_opf = solver.solve(model_multi_opf, tee=True) - - if check_status(results_multi_opf): - logging.info(f'Objective: {model_multi_opf.obj():.4e}') - - # %% Process results - results = nygrid_sim.get_results_multi_opf(model_multi_opf) - logging.info(f'Cost: {results["COST"]:.4e}') - - # Format thermal generation results - results_pg = results['PG'] - thermal_pg = results_pg.iloc[:, :num_thermal] - thermal_pg = thermal_pg.rename(columns=gen_rename) - thermal_pg.index.name = 'TimeStamp' - - # Save thermal generation to CSV - filename = f'thermal_w_solar_{start_datetime.strftime("%Y%m%d%H")}_{end_datetime.strftime("%Y%m%d%H")}.csv' - thermal_pg.to_csv(os.path.join(results_dir, 'w_solar', filename)) - logging.info(f'Saved thermal generation results in {filename}') - - # Save simulation results to pickle - filename = f'nygrid_sim_w_solar_{start_datetime.strftime("%Y%m%d%H")}_{end_datetime.strftime("%Y%m%d%H")}.pkl' - with open(os.path.join(results_dir, 'w_solar', filename), 'wb') as f: - pickle.dump([nygrid_sim, model_multi_opf, results], f) - logging.info(f'Saved simulation results in {filename}') - elapsed = time.time() - t - logging.info(f'Elapsed time: {elapsed:.2f} seconds') - logging.info('-----------------------------------------------------------------') - - else: - logging.error('No solution found.') - logging.info('-----------------------------------------------------------------') diff --git a/examples/ex_opf_w_offshore_wind.py b/examples/ex_opf_w_offshore_wind.py deleted file mode 100644 index 43c8bfb..0000000 --- a/examples/ex_opf_w_offshore_wind.py +++ /dev/null @@ -1,260 +0,0 @@ -""" -Run multi-period OPF with 2018 data -with renewable generators -including wind and solar - -""" -# %% Packages -import os -from pyomo.environ import * -import numpy as np -import pandas as pd -import matplotlib.pyplot as plt -from datetime import datetime, timedelta -from nygrid.nygrid import NYGrid -import pickle -import time - -t = time.time() - -# %% Set up directories -cwd = os.getcwd() -if 'examples' in cwd: - parent_dir = os.path.dirname(cwd) - data_dir = os.path.join(parent_dir, 'data') -else: - data_dir = os.path.join(cwd, 'data') - -grid_data_dir = os.path.join(data_dir, 'grid') -if not os.path.exists(grid_data_dir): - raise FileNotFoundError('Grid data directory not found.') - -print('Grid data directory: {}'.format(grid_data_dir)) - -fig_dir = os.path.join(os.path.dirname(data_dir), 'figures') -print('Figure directory: {}'.format(fig_dir)) - -results_dir = os.path.join(os.path.dirname(data_dir), 'results') -print('Results directory: {}'.format(results_dir)) - -if not os.path.exists(os.path.join(results_dir, 'w_offshore')): - os.mkdir(os.path.join(results_dir, 'w_offshore')) - -solar_data_dir = os.path.join(data_dir, 'solar') -print('Solar data directory: {}'.format(solar_data_dir)) - -onshore_wind_data_dir = os.path.join(data_dir, 'onshore_wind') -print('Onshore wind data directory: {}'.format(onshore_wind_data_dir)) - -offshore_wind_data_dir = os.path.join(data_dir, 'offshore_wind') -print('Offshore wind data directory: {}'.format(offshore_wind_data_dir)) - -# %% Read grid data -start_date = datetime(2018, 1, 1, 0, 0, 0) -end_date = datetime(2019, 1, 1, 0, 0, 0) -timestamp_list = pd.date_range(start_date, end_date, freq='1D') - -# Read load profile -load_profile = pd.read_csv(os.path.join(grid_data_dir, f'load_profile_{start_date.year}.csv'), - parse_dates=['TimeStamp'], index_col='TimeStamp') -load_profile.index.freq = 'H' - -# Read generation profile -gen_profile = pd.read_csv(os.path.join(grid_data_dir, f'gen_profile_{start_date.year}.csv'), - parse_dates=['TimeStamp'], index_col='TimeStamp') -gen_profile.index.freq = 'H' - -# Read generator capacity limit profile -genmax_profile = pd.read_csv(os.path.join(grid_data_dir, f'genmax_profile_{start_date.year}.csv'), - parse_dates=['TimeStamp'], index_col='TimeStamp') -genmax_profile.index.freq = 'H' - -genmin_profile = pd.read_csv(os.path.join(grid_data_dir, f'genmin_profile_{start_date.year}.csv'), - parse_dates=['TimeStamp'], index_col='TimeStamp') -genmin_profile.index.freq = 'H' - -# Read generator ramp rate profile -genramp30_profile = pd.read_csv(os.path.join(grid_data_dir, f'genramp30_profile_{start_date.year}.csv'), - parse_dates=['TimeStamp'], index_col='TimeStamp') -genramp30_profile.index.freq = 'H' - -# Read generator cost profile (linear) -gencost0_profile = pd.read_csv(os.path.join(grid_data_dir, f'gencost0_profile_{start_date.year}.csv'), - parse_dates=['TimeStamp'], index_col='TimeStamp') -gencost0_profile.index.freq = 'H' - -gencost1_profile = pd.read_csv(os.path.join(grid_data_dir, f'gencost1_profile_{start_date.year}.csv'), - parse_dates=['TimeStamp'], index_col='TimeStamp') -gencost1_profile.index.freq = 'H' - -# %% Read renewable generation data -# Renewable generation time series -current_solar_gen = pd.read_csv(os.path.join(solar_data_dir, f'current_solar_gen_1hr.csv'), - parse_dates=['Time'], index_col='Time') -current_solar_gen.index.freq = 'H' -current_solar_gen.columns = current_solar_gen.columns.astype(int) - -future_solar_gen = pd.read_csv(os.path.join(solar_data_dir, f'future_solar_gen_1hr.csv'), - parse_dates=['Time'], index_col='Time') -future_solar_gen.index.freq = 'H' -future_solar_gen.columns = future_solar_gen.columns.astype(int) - -onshore_wind_gen = pd.read_csv(os.path.join(onshore_wind_data_dir, f'current_wind_gen_1hr.csv'), - parse_dates=['Time'], index_col='Time') -onshore_wind_gen.index.freq = 'H' -onshore_wind_gen.columns = onshore_wind_gen.columns.astype(int) - -offshore_wind_gen = pd.read_csv(os.path.join(offshore_wind_data_dir, f'power_load_2018.csv'), - parse_dates=['timestamp'], index_col='timestamp') -offshore_wind_gen.index = offshore_wind_gen.index.tz_localize('US/Eastern', ambiguous='infer') -offshore_wind_gen.index.freq = 'H' - -# Renewable allocation table -current_solar_2bus = pd.read_csv(os.path.join(solar_data_dir, f'solar_farms_2bus.csv'), - index_col='zip_code') -future_solar_2bus = pd.read_csv(os.path.join(solar_data_dir, f'future_solar_farms_2bus.csv'), index_col=0) -onshore_wind_2bus = pd.read_csv(os.path.join(onshore_wind_data_dir, f'onshore_wind_2bus.csv'), index_col=0) - -# Aggregate current solar generation -groupby_dict = current_solar_2bus['busIdx'].to_dict() -current_solar_gen_agg = current_solar_gen.groupby(groupby_dict, axis=1).sum() -current_solar_gen_agg = current_solar_gen_agg/1e3 # convert from kW to MW - -# Aggregate future solar generation -groupby_dict = future_solar_2bus['busIdx'].to_dict() -future_solar_gen_agg = future_solar_gen.groupby(groupby_dict, axis=1).sum() -future_solar_gen_agg = future_solar_gen_agg/1e3 # convert from kW to MW - -# Aggregate onshore wind generation -groupby_dict = onshore_wind_2bus['busIdx'].to_dict() -onshore_wind_gen_agg = onshore_wind_gen.groupby(groupby_dict, axis=1).sum() -onshore_wind_gen_agg = onshore_wind_gen_agg/1e3 # convert from kW to MW - -# Aggregate offshore wind generation -offshore_wind_gen_agg = offshore_wind_gen[['power_nyc', 'power_li']].rename( - columns={'power_nyc': 81, 'power_li': 79}) - -# %% Update load profile -# Tread renewable generation as negative load -load_profile_copy = load_profile.copy() -load_profile_copy.columns = [int(col.replace('Bus', '')) for col in load_profile_copy.columns] -load_profile_copy.index = current_solar_gen_agg.index - -# 18.08% of current solar generation is built before 2018 (base year) -# Scale down current solar generation by 18.08% -pct_current_solar_built = 0.1808 -current_solar_gen_agg = current_solar_gen_agg * (1-pct_current_solar_built) - -# 90.6% of onshore wind generation is built before 2018 (base year) -# Scale down onshore wind generation by 90.6% -pct_onshore_wind_built = 0.906 -onshore_wind_gen_agg = onshore_wind_gen_agg * (1-pct_onshore_wind_built) - -# Total renewable generation by bus -total_renewable = pd.DataFrame() -total_renewable = total_renewable.add(onshore_wind_gen_agg, fill_value=0) -total_renewable = total_renewable.add(offshore_wind_gen_agg, fill_value=0) -# total_renewable = total_renewable.add(current_solar_gen_agg, fill_value=0) -# total_renewable = total_renewable.add(future_solar_gen_agg, fill_value=0) - -# Load profile after subtracting renewable generation -load_profile_renewable = load_profile_copy.subtract(total_renewable, fill_value=0) - -# Scale down external load by the same percentage as NYS -ny_buses = np.arange(37, 83) -ext_buses = np.array(list(set(load_profile_copy.columns) - set(ny_buses))) -ny_load_profile = load_profile_copy[ny_buses] -ext_load_profile = load_profile_copy[ext_buses] - -# NY load change rate in each hour -ny_change_pct = ((load_profile_renewable[ny_buses].sum(axis=1) - - ny_load_profile.sum(axis=1)) / ny_load_profile.sum(axis=1)) -print(f'NYS annual load change: {ny_change_pct.mean()*100:.2f}%') - -ext_change = ext_load_profile.multiply(ny_change_pct*0.2, axis=0) -load_profile_renewable = load_profile_renewable.add(ext_change, fill_value=0) -print(f'External load change is set to 20% of NYS load change.') - -# Reset timestamp index to time zone unaware -load_profile_copy.index = load_profile.index -load_profile_renewable.index = load_profile.index - -# %% Read thermal generator info table -filename = os.path.join(data_dir, 'genInfo.csv') -gen_info = pd.read_csv(filename) -num_thermal = gen_info.shape[0] -gen_rename = {gen_info.index[i]: gen_info.NYISOName[i] for i in range(num_thermal)} - -# %% Set up OPF model -timestamp_list = pd.date_range(start_date, end_date, freq='1D') - -# Loop through all days -for d in range(len(timestamp_list)-1): - # Run one day at a time - start_datetime = timestamp_list[d] - end_datetime = start_datetime + timedelta(hours=23) - - # Read MATPOWER case file - ppc_filename = os.path.join(data_dir, 'ny_grid.mat') - - nygrid_sim = NYGrid(ppc_filename, - start_datetime=start_datetime.strftime('%m-%d-%Y %H'), - end_datetime=end_datetime.strftime('%m-%d-%Y %H'), - verbose=True) - - # Read grid data - nygrid_sim.set_load_sch(load_profile_renewable) - nygrid_sim.set_gen_mw_sch(gen_profile) - nygrid_sim.set_gen_max_sch(genmax_profile) - nygrid_sim.set_gen_min_sch(genmin_profile) - nygrid_sim.set_gen_ramp_sch(genramp30_profile) - nygrid_sim.set_gen_cost_sch(gencost0_profile, gencost1_profile) - - # Process ppc - nygrid_sim.process_ppc() - - # Set generator initial condition - if d == 0: - last_gen = None - print('No initial condition.') - else: - last_gen = nygrid_sim.get_last_gen(model_multi_opf) - print('Initial condition set from previous day.') - - nygrid_sim.set_gen_init_data(gen_init=last_gen) - - # Check input - nygrid_sim.check_input_dim() - - # Initialize single period OPF - model_multi_opf = nygrid_sim.create_multi_opf_soft(slack_cost_weight=1e21) - - solver = SolverFactory('gurobi') - results_multi_opf = solver.solve(model_multi_opf, tee=True) - - if check_status(results_multi_opf): - print(f'Objective: {model_multi_opf.obj():.4e}') - - # %% Process results - results = nygrid_sim.get_results_multi_opf(model_multi_opf) - print(f'Cost: {results["COST"]:.4e}') - - # Format thermal generation results - results_pg = results['PG'] - thermal_pg = results_pg.iloc[:, :num_thermal] - thermal_pg = thermal_pg.rename(columns=gen_rename) - thermal_pg.index.name = 'TimeStamp' - - # Save thermal generation to CSV - filename = f'thermal_w_offshore_{start_datetime.strftime("%Y%m%d%H")}_{end_datetime.strftime("%Y%m%d%H")}.csv' - thermal_pg.to_csv(os.path.join(results_dir, 'w_offshore', filename)) - print(f'Saved thermal generation results in {filename}') - - # Save simulation results to pickle - filename = f'nygrid_sim_w_offshore_{start_datetime.strftime("%Y%m%d%H")}_{end_datetime.strftime("%Y%m%d%H")}.pkl' - with open(os.path.join(results_dir, 'w_offshore', filename), 'wb') as f: - pickle.dump([nygrid_sim, model_multi_opf, results], f) - print(f'Saved simulation results in {filename}') - elapsed = time.time() - t - print(f'Elapsed time: {elapsed:.2f} seconds') - print('-----------------------------------------------------------------') diff --git a/examples/ex_opf_w_renew.py b/examples/ex_opf_w_renew.py deleted file mode 100644 index f4b6aab..0000000 --- a/examples/ex_opf_w_renew.py +++ /dev/null @@ -1,260 +0,0 @@ -""" -Run multi-period OPF with 2018 data -with renewable generators -including wind and solar - -""" -# %% Packages -import os -from pyomo.environ import * -import numpy as np -import pandas as pd -import matplotlib.pyplot as plt -from datetime import datetime, timedelta -from nygrid.nygrid import NYGrid -import pickle -import time - -t = time.time() - -# %% Set up directories -cwd = os.getcwd() -if 'examples' in cwd: - parent_dir = os.path.dirname(cwd) - data_dir = os.path.join(parent_dir, 'data') -else: - data_dir = os.path.join(cwd, 'data') - -grid_data_dir = os.path.join(data_dir, 'grid') -if not os.path.exists(grid_data_dir): - raise FileNotFoundError('Grid data directory not found.') - -print('Grid data directory: {}'.format(grid_data_dir)) - -fig_dir = os.path.join(os.path.dirname(data_dir), 'figures') -print('Figure directory: {}'.format(fig_dir)) - -results_dir = os.path.join(os.path.dirname(data_dir), 'results') -print('Results directory: {}'.format(results_dir)) - -if not os.path.exists(os.path.join(results_dir, 'w_renew')): - os.mkdir(os.path.join(results_dir, 'w_renew')) - -solar_data_dir = os.path.join(data_dir, 'solar') -print('Solar data directory: {}'.format(solar_data_dir)) - -onshore_wind_data_dir = os.path.join(data_dir, 'onshore_wind') -print('Onshore wind data directory: {}'.format(onshore_wind_data_dir)) - -offshore_wind_data_dir = os.path.join(data_dir, 'offshore_wind') -print('Offshore wind data directory: {}'.format(offshore_wind_data_dir)) - -# %% Read grid data -start_date = datetime(2018, 1, 1, 0, 0, 0) -end_date = datetime(2019, 1, 1, 0, 0, 0) -timestamp_list = pd.date_range(start_date, end_date, freq='1D') - -# Read load profile -load_profile = pd.read_csv(os.path.join(grid_data_dir, f'load_profile_{start_date.year}.csv'), - parse_dates=['TimeStamp'], index_col='TimeStamp') -load_profile.index.freq = 'H' - -# Read generation profile -gen_profile = pd.read_csv(os.path.join(grid_data_dir, f'gen_profile_{start_date.year}.csv'), - parse_dates=['TimeStamp'], index_col='TimeStamp') -gen_profile.index.freq = 'H' - -# Read generator capacity limit profile -genmax_profile = pd.read_csv(os.path.join(grid_data_dir, f'genmax_profile_{start_date.year}.csv'), - parse_dates=['TimeStamp'], index_col='TimeStamp') -genmax_profile.index.freq = 'H' - -genmin_profile = pd.read_csv(os.path.join(grid_data_dir, f'genmin_profile_{start_date.year}.csv'), - parse_dates=['TimeStamp'], index_col='TimeStamp') -genmin_profile.index.freq = 'H' - -# Read generator ramp rate profile -genramp30_profile = pd.read_csv(os.path.join(grid_data_dir, f'genramp30_profile_{start_date.year}.csv'), - parse_dates=['TimeStamp'], index_col='TimeStamp') -genramp30_profile.index.freq = 'H' - -# Read generator cost profile (linear) -gencost0_profile = pd.read_csv(os.path.join(grid_data_dir, f'gencost0_profile_{start_date.year}.csv'), - parse_dates=['TimeStamp'], index_col='TimeStamp') -gencost0_profile.index.freq = 'H' - -gencost1_profile = pd.read_csv(os.path.join(grid_data_dir, f'gencost1_profile_{start_date.year}.csv'), - parse_dates=['TimeStamp'], index_col='TimeStamp') -gencost1_profile.index.freq = 'H' - -# %% Read renewable generation data -# Renewable generation time series -current_solar_gen = pd.read_csv(os.path.join(solar_data_dir, f'current_solar_gen_1hr.csv'), - parse_dates=['Time'], index_col='Time') -current_solar_gen.index.freq = 'H' -current_solar_gen.columns = current_solar_gen.columns.astype(int) - -future_solar_gen = pd.read_csv(os.path.join(solar_data_dir, f'future_solar_gen_1hr.csv'), - parse_dates=['Time'], index_col='Time') -future_solar_gen.index.freq = 'H' -future_solar_gen.columns = future_solar_gen.columns.astype(int) - -onshore_wind_gen = pd.read_csv(os.path.join(onshore_wind_data_dir, f'current_wind_gen_1hr.csv'), - parse_dates=['Time'], index_col='Time') -onshore_wind_gen.index.freq = 'H' -onshore_wind_gen.columns = onshore_wind_gen.columns.astype(int) - -offshore_wind_gen = pd.read_csv(os.path.join(offshore_wind_data_dir, f'power_load_2018.csv'), - parse_dates=['timestamp'], index_col='timestamp') -offshore_wind_gen.index = offshore_wind_gen.index.tz_localize('US/Eastern', ambiguous='infer') -offshore_wind_gen.index.freq = 'H' - -# Renewable allocation table -current_solar_2bus = pd.read_csv(os.path.join(solar_data_dir, f'solar_farms_2bus.csv'), - index_col='zip_code') -future_solar_2bus = pd.read_csv(os.path.join(solar_data_dir, f'future_solar_farms_2bus.csv'), index_col=0) -onshore_wind_2bus = pd.read_csv(os.path.join(onshore_wind_data_dir, f'onshore_wind_2bus.csv'), index_col=0) - -# Aggregate current solar generation -groupby_dict = current_solar_2bus['busIdx'].to_dict() -current_solar_gen_agg = current_solar_gen.groupby(groupby_dict, axis=1).sum() -current_solar_gen_agg = current_solar_gen_agg/1e3 # convert from kW to MW - -# Aggregate future solar generation -groupby_dict = future_solar_2bus['busIdx'].to_dict() -future_solar_gen_agg = future_solar_gen.groupby(groupby_dict, axis=1).sum() -future_solar_gen_agg = future_solar_gen_agg/1e3 # convert from kW to MW - -# Aggregate onshore wind generation -groupby_dict = onshore_wind_2bus['busIdx'].to_dict() -onshore_wind_gen_agg = onshore_wind_gen.groupby(groupby_dict, axis=1).sum() -onshore_wind_gen_agg = onshore_wind_gen_agg/1e3 # convert from kW to MW - -# Aggregate offshore wind generation -offshore_wind_gen_agg = offshore_wind_gen[['power_nyc', 'power_li']].rename( - columns={'power_nyc': 81, 'power_li': 79}) - -# %% Update load profile -# Tread renewable generation as negative load -load_profile_copy = load_profile.copy() -load_profile_copy.columns = [int(col.replace('Bus', '')) for col in load_profile_copy.columns] -load_profile_copy.index = current_solar_gen_agg.index - -# 18.08% of current solar generation is built before 2018 (base year) -# Scale down current solar generation by 18.08% -pct_current_solar_built = 0.1808 -current_solar_gen_agg = current_solar_gen_agg * (1-pct_current_solar_built) - -# 90.6% of onshore wind generation is built before 2018 (base year) -# Scale down onshore wind generation by 90.6% -pct_onshore_wind_built = 0.906 -onshore_wind_gen_agg = onshore_wind_gen_agg * (1-pct_onshore_wind_built) - -# Total renewable generation by bus -total_renewable = pd.DataFrame() -total_renewable = total_renewable.add(onshore_wind_gen_agg, fill_value=0) -total_renewable = total_renewable.add(offshore_wind_gen_agg, fill_value=0) -total_renewable = total_renewable.add(current_solar_gen_agg, fill_value=0) -total_renewable = total_renewable.add(future_solar_gen_agg, fill_value=0) - -# Load profile after subtracting renewable generation -load_profile_renewable = load_profile_copy.subtract(total_renewable, fill_value=0) - -# Scale down external load by the same percentage as NYS -ny_buses = np.arange(37, 83) -ext_buses = np.array(list(set(load_profile_copy.columns) - set(ny_buses))) -ny_load_profile = load_profile_copy[ny_buses] -ext_load_profile = load_profile_copy[ext_buses] - -# NY load change rate in each hour -ny_change_pct = ((load_profile_renewable[ny_buses].sum(axis=1) - - ny_load_profile.sum(axis=1)) / ny_load_profile.sum(axis=1)) -print(f'NYS annual load change: {ny_change_pct.mean()*100:.2f}%') - -ext_change = ext_load_profile.multiply(ny_change_pct*0.2, axis=0) -load_profile_renewable = load_profile_renewable.add(ext_change, fill_value=0) -print(f'External load change is set to 20% of NYS load change.') - -# Reset timestamp index to time zone unaware -load_profile_copy.index = load_profile.index -load_profile_renewable.index = load_profile.index - -# %% Read thermal generator info table -filename = os.path.join(data_dir, 'genInfo.csv') -gen_info = pd.read_csv(filename) -num_thermal = gen_info.shape[0] -gen_rename = {gen_info.index[i]: gen_info.NYISOName[i] for i in range(num_thermal)} - -# %% Set up OPF model -timestamp_list = pd.date_range(start_date, end_date, freq='1D') - -# Loop through all days -for d in range(len(timestamp_list)-1): - # Run one day at a time - start_datetime = timestamp_list[d] - end_datetime = start_datetime + timedelta(hours=23) - - # Read MATPOWER case file - ppc_filename = os.path.join(data_dir, 'ny_grid.mat') - - nygrid_sim = NYGrid(ppc_filename, - start_datetime=start_datetime.strftime('%m-%d-%Y %H'), - end_datetime=end_datetime.strftime('%m-%d-%Y %H'), - verbose=True) - - # Read grid data - nygrid_sim.set_load_sch(load_profile_renewable) - nygrid_sim.set_gen_mw_sch(gen_profile) - nygrid_sim.set_gen_max_sch(genmax_profile) - nygrid_sim.set_gen_min_sch(genmin_profile) - nygrid_sim.set_gen_ramp_sch(genramp30_profile) - nygrid_sim.set_gen_cost_sch(gencost0_profile, gencost1_profile) - - # Process ppc - nygrid_sim.process_ppc() - - # Set generator initial condition - if d == 0: - last_gen = None - print('No initial condition.') - else: - last_gen = nygrid_sim.get_last_gen(model_multi_opf) - print('Initial condition set from previous day.') - - nygrid_sim.set_gen_init_data(gen_init=last_gen) - - # Check input - nygrid_sim.check_input_dim() - - # Initialize single period OPF - model_multi_opf = nygrid_sim.create_multi_opf_soft(slack_cost_weight=1e21) - - solver = SolverFactory('gurobi') - results_multi_opf = solver.solve(model_multi_opf, tee=True) - - if check_status(results_multi_opf): - print(f'Objective: {model_multi_opf.obj():.4e}') - - # %% Process results - results = nygrid_sim.get_results_multi_opf(model_multi_opf) - print(f'Cost: {results["COST"]:.4e}') - - # Format thermal generation results - results_pg = results['PG'] - thermal_pg = results_pg.iloc[:, :num_thermal] - thermal_pg = thermal_pg.rename(columns=gen_rename) - thermal_pg.index.name = 'TimeStamp' - - # Save thermal generation to CSV - filename = f'thermal_w_renew_{start_datetime.strftime("%Y%m%d%H")}_{end_datetime.strftime("%Y%m%d%H")}.csv' - thermal_pg.to_csv(os.path.join(results_dir, 'w_renew', filename)) - print(f'Saved thermal generation results in {filename}') - - # Save simulation results to pickle - filename = f'nygrid_sim_w_renew_{start_datetime.strftime("%Y%m%d%H")}_{end_datetime.strftime("%Y%m%d%H")}.pkl' - with open(os.path.join(results_dir, 'w_renew', filename), 'wb') as f: - pickle.dump([nygrid_sim, model_multi_opf, results], f) - print(f'Saved simulation results in {filename}') - elapsed = time.time() - t - print(f'Elapsed time: {elapsed:.2f} seconds') - print('-----------------------------------------------------------------') diff --git a/examples/ex_opf_w_renew_elec.py b/examples/ex_opf_w_renew_elec.py deleted file mode 100644 index 702ec46..0000000 --- a/examples/ex_opf_w_renew_elec.py +++ /dev/null @@ -1,310 +0,0 @@ -""" -Run multi-period OPF with 2018 data -with renewable generators -including wind and solar - -""" -# %% Packages -import os -from pyomo.environ import * -import numpy as np -import pandas as pd -import matplotlib.pyplot as plt -from datetime import datetime, timedelta -from nygrid.nygrid import NYGrid -import pickle -import time - -t = time.time() - -# %% Set up directories -cwd = os.getcwd() -if 'examples' in cwd: - parent_dir = os.path.dirname(cwd) - data_dir = os.path.join(parent_dir, 'data') -else: - data_dir = os.path.join(cwd, 'data') - -grid_data_dir = os.path.join(data_dir, 'grid') -if not os.path.exists(grid_data_dir): - raise FileNotFoundError('Grid data directory not found.') - -print('Grid data directory: {}'.format(grid_data_dir)) - -fig_dir = os.path.join(os.path.dirname(data_dir), 'figures') -print('Figure directory: {}'.format(fig_dir)) - -results_dir = os.path.join(os.path.dirname(data_dir), 'results') -print('Results directory: {}'.format(results_dir)) - -if not os.path.exists(os.path.join(results_dir, 'w_renew_elec')): - os.mkdir(os.path.join(results_dir, 'w_renew_elec')) - -solar_data_dir = os.path.join(data_dir, 'solar') -print('Solar data directory: {}'.format(solar_data_dir)) - -onshore_wind_data_dir = os.path.join(data_dir, 'onshore_wind') -print('Onshore wind data directory: {}'.format(onshore_wind_data_dir)) - -offshore_wind_data_dir = os.path.join(data_dir, 'offshore_wind') -print('Offshore wind data directory: {}'.format(offshore_wind_data_dir)) - -buildings_data_dir = os.path.join(data_dir, 'buildings') -print('Buildings data directory: {}'.format(buildings_data_dir)) - -# %% Read grid data -start_date = datetime(2018, 1, 1, 0, 0, 0) -end_date = datetime(2019, 1, 1, 0, 0, 0) -timestamp_list = pd.date_range(start_date, end_date, freq='1D') - -# Read load profile -load_profile = pd.read_csv(os.path.join(grid_data_dir, f'load_profile_{start_date.year}.csv'), - parse_dates=['TimeStamp'], index_col='TimeStamp') -load_profile.index.freq = 'H' - -# Read generation profile -gen_profile = pd.read_csv(os.path.join(grid_data_dir, f'gen_profile_{start_date.year}.csv'), - parse_dates=['TimeStamp'], index_col='TimeStamp') -gen_profile.index.freq = 'H' - -# Read generator capacity limit profile -genmax_profile = pd.read_csv(os.path.join(grid_data_dir, f'genmax_profile_{start_date.year}.csv'), - parse_dates=['TimeStamp'], index_col='TimeStamp') -genmax_profile.index.freq = 'H' - -genmin_profile = pd.read_csv(os.path.join(grid_data_dir, f'genmin_profile_{start_date.year}.csv'), - parse_dates=['TimeStamp'], index_col='TimeStamp') -genmin_profile.index.freq = 'H' - -# Read generator ramp rate profile -genramp30_profile = pd.read_csv(os.path.join(grid_data_dir, f'genramp30_profile_{start_date.year}.csv'), - parse_dates=['TimeStamp'], index_col='TimeStamp') -genramp30_profile.index.freq = 'H' - -# Read generator cost profile (linear) -gencost0_profile = pd.read_csv(os.path.join(grid_data_dir, f'gencost0_profile_{start_date.year}.csv'), - parse_dates=['TimeStamp'], index_col='TimeStamp') -gencost0_profile.index.freq = 'H' - -gencost1_profile = pd.read_csv(os.path.join(grid_data_dir, f'gencost1_profile_{start_date.year}.csv'), - parse_dates=['TimeStamp'], index_col='TimeStamp') -gencost1_profile.index.freq = 'H' - -# %% Read renewable generation data -# Renewable generation time series -current_solar_gen = pd.read_csv(os.path.join(solar_data_dir, f'current_solar_gen_1hr.csv'), - parse_dates=['Time'], index_col='Time') -current_solar_gen.index.freq = 'H' -current_solar_gen.columns = current_solar_gen.columns.astype(int) - -future_solar_gen = pd.read_csv(os.path.join(solar_data_dir, f'future_solar_gen_1hr.csv'), - parse_dates=['Time'], index_col='Time') -future_solar_gen.index.freq = 'H' -future_solar_gen.columns = future_solar_gen.columns.astype(int) - -onshore_wind_gen = pd.read_csv(os.path.join(onshore_wind_data_dir, f'current_wind_gen_1hr.csv'), - parse_dates=['Time'], index_col='Time') -onshore_wind_gen.index.freq = 'H' -onshore_wind_gen.columns = onshore_wind_gen.columns.astype(int) - -offshore_wind_gen = pd.read_csv(os.path.join(offshore_wind_data_dir, f'power_load_2018.csv'), - parse_dates=['timestamp'], index_col='timestamp') -offshore_wind_gen.index = offshore_wind_gen.index.tz_localize('US/Eastern', ambiguous='infer') -offshore_wind_gen.index.freq = 'H' - -# Renewable allocation table -current_solar_2bus = pd.read_csv(os.path.join(solar_data_dir, f'solar_farms_2bus.csv'), - index_col='zip_code') -future_solar_2bus = pd.read_csv(os.path.join(solar_data_dir, f'future_solar_farms_2bus.csv'), index_col=0) -onshore_wind_2bus = pd.read_csv(os.path.join(onshore_wind_data_dir, f'onshore_wind_2bus.csv'), index_col=0) - -# Aggregate current solar generation -groupby_dict = current_solar_2bus['busIdx'].to_dict() -current_solar_gen_agg = current_solar_gen.groupby(groupby_dict, axis=1).sum() -current_solar_gen_agg = current_solar_gen_agg/1e3 # convert from kW to MW - -# Aggregate future solar generation -groupby_dict = future_solar_2bus['busIdx'].to_dict() -future_solar_gen_agg = future_solar_gen.groupby(groupby_dict, axis=1).sum() -future_solar_gen_agg = future_solar_gen_agg/1e3 # convert from kW to MW - -# Aggregate onshore wind generation -groupby_dict = onshore_wind_2bus['busIdx'].to_dict() -onshore_wind_gen_agg = onshore_wind_gen.groupby(groupby_dict, axis=1).sum() -onshore_wind_gen_agg = onshore_wind_gen_agg/1e3 # convert from kW to MW - -# Aggregate offshore wind generation -offshore_wind_gen_agg = offshore_wind_gen[['power_nyc', 'power_li']].rename( - columns={'power_nyc': 81, 'power_li': 79}) - -# %% Read electrification data -# Read county to bus allocation table -county_2bus = pd.read_csv(os.path.join(buildings_data_dir, 'county_centroids_2bus.csv'), - index_col=0) -fips_list = ['G'+s[:2]+'0'+s[2:]+'0' for s in county_2bus['FIPS_CODE'].astype(str)] - -# County-level load data -# Residential buildings current load (upgrade=0) -pickle_dir = os.path.join(buildings_data_dir, 'euss_processed', 'upgrade=0') -res_current_load = pd.DataFrame() - -for ii in range(len(fips_list)): - pickle_file = os.path.join(pickle_dir, f'{fips_list[0]}_elec_total.pkl') - # Read pickle file - with open(os.path.join(pickle_dir, pickle_file), 'rb') as f: - county_load = pickle.load(f) - county_load = county_load.rename(columns={'total': fips_list[ii]}) - res_current_load = pd.concat([res_current_load, county_load], axis=1) - -# Residential buildings future load (upgrade=10) -pickle_dir = os.path.join(buildings_data_dir, 'euss_processed', 'upgrade=10') -res_elec_load = pd.DataFrame() - -for ii in range(len(fips_list)): - pickle_file = os.path.join(pickle_dir, f'{fips_list[0]}_elec_total.pkl') - # Read pickle file - with open(os.path.join(pickle_dir, pickle_file), 'rb') as f: - county_load = pickle.load(f) - county_load = county_load.rename(columns={'total': fips_list[ii]}) - res_elec_load = pd.concat([res_elec_load, county_load], axis=1) - -# Changes in residential load -res_load_change = res_elec_load - res_current_load -res_load_change.index.name = 'Time' - -# Aggregate building load to buses -groupby_dict = dict(zip(fips_list, county_2bus['busIdx'])) -res_load_change_agg = res_load_change.groupby(groupby_dict, axis=1).sum() -res_load_change_agg = res_load_change_agg/1e3 # convert from kW to MW - -# Scale down residential load change by 50% -res_load_change_agg = res_load_change_agg*0.5 - -# %% Update load profile -# Tread renewable generation as negative load -load_profile_copy = load_profile.copy() -load_profile_copy.columns = [int(col.replace('Bus', '')) for col in load_profile_copy.columns] -load_profile_copy.index = current_solar_gen_agg.index - -# 18.08% of current solar generation is built before 2018 (base year) -# Scale down current solar generation by 18.08% -pct_current_solar_built = 0.1808 -current_solar_gen_agg = current_solar_gen_agg * (1-pct_current_solar_built) - -# 90.6% of onshore wind generation is built before 2018 (base year) -# Scale down onshore wind generation by 90.6% -pct_onshore_wind_built = 0.906 -onshore_wind_gen_agg = onshore_wind_gen_agg * (1-pct_onshore_wind_built) - -# Total renewable generation by bus -total_renewable = pd.DataFrame() -total_renewable = total_renewable.add(onshore_wind_gen_agg, fill_value=0) -total_renewable = total_renewable.add(offshore_wind_gen_agg, fill_value=0) -total_renewable = total_renewable.add(current_solar_gen_agg, fill_value=0) -total_renewable = total_renewable.add(future_solar_gen_agg, fill_value=0) - -# Load profile after subtracting renewable generation -load_profile_renewable = load_profile_copy.subtract(total_renewable, fill_value=0) - -# Load profile with residential load change -res_load_change_agg.index = current_solar_gen_agg.index -load_profile_renewable = load_profile_renewable.add(res_load_change_agg, fill_value=0) - -# Scale down external load by the same percentage as NYS -ny_buses = np.arange(37, 83) -ext_buses = np.array(list(set(load_profile_copy.columns) - set(ny_buses))) -ny_load_profile = load_profile_copy[ny_buses] -ext_load_profile = load_profile_copy[ext_buses] - -# NY load change rate in each hour -ny_change_pct = ((load_profile_renewable[ny_buses].sum(axis=1) - - ny_load_profile.sum(axis=1)) / ny_load_profile.sum(axis=1)) -print(f'NYS annual load change: {ny_change_pct.mean()*100:.2f}%') - -ext_change = ext_load_profile.multiply(ny_change_pct*0.2, axis=0) -load_profile_renewable = load_profile_renewable.add(ext_change, fill_value=0) -print(f'External load change is set to 20% of NYS load change.') - -# Reset timestamp index to time zone unaware -load_profile_copy.index = load_profile.index -load_profile_renewable.index = load_profile.index - -# %% Read thermal generator info table -filename = os.path.join(data_dir, 'genInfo.csv') -gen_info = pd.read_csv(filename) -num_thermal = gen_info.shape[0] -gen_rename = {gen_info.index[i]: gen_info.NYISOName[i] for i in range(num_thermal)} - -# %% Set up OPF model -timestamp_list = pd.date_range(start_date, end_date, freq='1D') - -# Loop through all days -for d in range(len(timestamp_list)-1): - # Run one day at a time - start_datetime = timestamp_list[d] - end_datetime = start_datetime + timedelta(hours=23) - - # Read MATPOWER case file - ppc_filename = os.path.join(data_dir, 'ny_grid.mat') - - nygrid_sim = NYGrid(ppc_filename, - start_datetime=start_datetime.strftime('%m-%d-%Y %H'), - end_datetime=end_datetime.strftime('%m-%d-%Y %H'), - verbose=True) - - # Read grid data - nygrid_sim.set_load_sch(load_profile_renewable) - nygrid_sim.set_gen_mw_sch(gen_profile) - nygrid_sim.set_gen_max_sch(genmax_profile) - nygrid_sim.set_gen_min_sch(genmin_profile) - nygrid_sim.set_gen_ramp_sch(genramp30_profile) - nygrid_sim.set_gen_cost_sch(gencost0_profile, gencost1_profile) - - # Process ppc - nygrid_sim.process_ppc() - - # Set generator initial condition - if d == 0: - last_gen = None - print('No initial condition.') - else: - last_gen = nygrid_sim.get_last_gen(model_multi_opf) - print('Initial condition set from previous day.') - - nygrid_sim.set_gen_init_data(gen_init=last_gen) - - # Check input - nygrid_sim.check_input_dim() - - # Initialize single period OPF - model_multi_opf = nygrid_sim.create_multi_opf_soft(slack_cost_weight=1e21) - - solver = SolverFactory('gurobi') - results_multi_opf = solver.solve(model_multi_opf, tee=True) - - if check_status(results_multi_opf): - print(f'Objective: {model_multi_opf.obj():.4e}') - - # %% Process results - results = nygrid_sim.get_results_multi_opf(model_multi_opf) - print(f'Cost: {results["COST"]:.4e}') - - # Format thermal generation results - results_pg = results['PG'] - thermal_pg = results_pg.iloc[:, :num_thermal] - thermal_pg = thermal_pg.rename(columns=gen_rename) - thermal_pg.index.name = 'TimeStamp' - - # Save thermal generation to CSV - filename = f'thermal_w_renew_elec_{start_datetime.strftime("%Y%m%d%H")}_{end_datetime.strftime("%Y%m%d%H")}.csv' - thermal_pg.to_csv(os.path.join(results_dir, 'w_renew_elec', filename)) - print(f'Saved thermal generation results in {filename}') - - # Save simulation results to pickle - filename = f'nygrid_sim_w_renew_elec_{start_datetime.strftime("%Y%m%d%H")}_{end_datetime.strftime("%Y%m%d%H")}.pkl' - with open(os.path.join(results_dir, 'w_renew_elec', filename), 'wb') as f: - pickle.dump([nygrid_sim, model_multi_opf, results], f) - print(f'Saved simulation results in {filename}') - elapsed = time.time() - t - print(f'Elapsed time: {elapsed:.2f} seconds') - print('-----------------------------------------------------------------') diff --git a/nygrid/run_nygrid.py b/nygrid/run_nygrid.py index 49baa76..f197bd8 100644 --- a/nygrid/run_nygrid.py +++ b/nygrid/run_nygrid.py @@ -1,7 +1,7 @@ import os - +import numpy as np import pandas as pd - +import pickle from nygrid.nygrid import NYGrid @@ -64,6 +64,168 @@ def read_grid_data(data_dir, year): return grid_data +def read_vre_data(solar_data_dir, onshore_wind_data_dir, offshore_wind_data_dir): + """ + + Parameters + ---------- + solar_data_dir: str + Directory of solar data + onshore_wind_data_dir: str + Directory of onshore wind data + offshore_wind_data_dir: str + Directory of offshore wind data + + Returns + ------- + vre_prop: pandas.DataFrame + VRE properties + genmax_profile_vre: pandas.DataFrame + VRE generation profiles + """ + + # Renewable generation time series + current_solar_gen = pd.read_csv(os.path.join(solar_data_dir, f'current_solar_gen_1hr.csv'), + parse_dates=['Time'], index_col='Time').asfreq('H') + current_solar_gen.columns = current_solar_gen.columns.astype(int) + + future_solar_gen = pd.read_csv(os.path.join(solar_data_dir, f'future_solar_gen_1hr.csv'), + parse_dates=['Time'], index_col='Time').asfreq('H') + future_solar_gen.columns = future_solar_gen.columns.astype(int) + + onshore_wind_gen = pd.read_csv(os.path.join(onshore_wind_data_dir, f'current_wind_gen_1hr.csv'), + parse_dates=['Time'], index_col='Time').asfreq('H') + onshore_wind_gen.columns = onshore_wind_gen.columns.astype(int) + + offshore_wind_gen = pd.read_csv(os.path.join(offshore_wind_data_dir, f'power_load_2018.csv'), + parse_dates=['timestamp'], index_col='timestamp') + offshore_wind_gen.index = offshore_wind_gen.index.tz_localize('US/Eastern', ambiguous='infer') + offshore_wind_gen.index.freq = 'H' + + # Wind farm capacity info + capacity = [816, 1260, 924, 1230] + capacity_nyc, capacity_li = np.sum(capacity[:2]), np.sum(capacity[2:]) + + # Correct offshore wind generation + offshore_wind_gen['power_nyc'] = np.where(offshore_wind_gen['power_nyc'] > capacity_nyc, capacity_nyc, + offshore_wind_gen['power_nyc']) + offshore_wind_gen['power_li'] = np.where(offshore_wind_gen['power_li'] > capacity_li, capacity_li, + offshore_wind_gen['power_li']) + + # Renewable allocation table + current_solar_2bus = pd.read_csv(os.path.join(solar_data_dir, f'solar_farms_2bus.csv'), index_col='zip_code') + future_solar_2bus = pd.read_csv(os.path.join(solar_data_dir, f'future_solar_farms_2bus.csv'), index_col=0) + onshore_wind_2bus = pd.read_csv(os.path.join(onshore_wind_data_dir, f'onshore_wind_2bus.csv'), index_col=0) + + # Aggregate current solar generation + groupby_dict = current_solar_2bus['busIdx'].to_dict() + current_solar_gen_agg = current_solar_gen.groupby(groupby_dict, axis=1).sum() + current_solar_gen_agg = current_solar_gen_agg / 1e3 # convert from kW to MW + + # Aggregate future solar generation + groupby_dict = future_solar_2bus['busIdx'].to_dict() + future_solar_gen_agg = future_solar_gen.groupby(groupby_dict, axis=1).sum() + future_solar_gen_agg = future_solar_gen_agg / 1e3 # convert from kW to MW + + # Aggregate onshore wind generation + groupby_dict = onshore_wind_2bus['busIdx'].to_dict() + onshore_wind_gen_agg = onshore_wind_gen.groupby(groupby_dict, axis=1).sum() + onshore_wind_gen_agg = onshore_wind_gen_agg / 1e3 # convert from kW to MW + + # Aggregate offshore wind generation + offshore_wind_gen_agg = offshore_wind_gen[['power_nyc', 'power_li']].rename( + columns={'power_nyc': 81, 'power_li': 82}) + + # 18.08% of current solar generation is built before 2018 (base year) + # Scale down current solar generation by 18.08% + pct_current_solar_built = 0.1808 + current_solar_gen_agg = current_solar_gen_agg * (1 - pct_current_solar_built) + + # 90.6% of onshore wind generation is built before 2018 (base year) + # Scale down onshore wind generation by 90.6% + pct_onshore_wind_built = 0.906 + onshore_wind_gen_agg = onshore_wind_gen_agg * (1 - pct_onshore_wind_built) + + vre_profiles = {'CurSol': current_solar_gen_agg, + 'FutSol': future_solar_gen_agg, + 'OnWind': onshore_wind_gen_agg, + 'OffWind': offshore_wind_gen_agg} + vre_prop_list = list() + for key, profile in vre_profiles.items(): + vre_prop_a = pd.DataFrame(data={'VRE_BUS': profile.columns, + 'VRE_PMAX': profile.max(axis=0), + 'VRE_PMIN': 0, + 'VRE_TYPE': 'Solar', + 'VRE_NAME': [f'{key}_{col}' for col in profile.columns]}) + vre_prop_list.append(vre_prop_a) + + # Combine gen_prop tables + vre_prop = pd.concat(vre_prop_list, ignore_index=True) + + # Combine genmax tables + genmax_profile_vre = pd.concat([current_solar_gen_agg, future_solar_gen_agg, + onshore_wind_gen_agg, offshore_wind_gen_agg], axis=1) + genmax_profile_vre.columns = vre_prop['VRE_NAME'] + genmax_profile_vre.index = genmax_profile_vre.index.tz_convert('US/Eastern').tz_localize(None) + + return vre_prop, genmax_profile_vre + + +def read_electrification_data(buildings_data_dir): + """ + + Parameters + ---------- + buildings_data_dir: str + Directory of buildings data + + Returns + ------- + res_load_change_bus: pandas.DataFrame + Changes in residential load + """ + # Read county to bus allocation table + county_2bus = pd.read_csv(os.path.join(buildings_data_dir, 'county_centroids_2bus.csv'), + index_col=0) + fips_list = ['G' + s[:2] + '0' + s[2:] + '0' for s in county_2bus['FIPS_CODE'].astype(str)] + + # County-level load data + # Residential buildings current load (upgrade=0) + pickle_dir = os.path.join(buildings_data_dir, 'euss_processed', 'upgrade=0') + res_current_load = pd.DataFrame() + + for ii in range(len(fips_list)): + pickle_file = os.path.join(pickle_dir, f'{fips_list[0]}_elec_total.pkl') + # Read pickle file + with open(os.path.join(pickle_dir, pickle_file), 'rb') as f: + county_load = pickle.load(f) + county_load = county_load.rename(columns={'total': fips_list[ii]}) + res_current_load = pd.concat([res_current_load, county_load], axis=1) + + # Residential buildings future load (upgrade=10) + pickle_dir = os.path.join(buildings_data_dir, 'euss_processed', 'upgrade=10') + res_elec_load = pd.DataFrame() + + for ii in range(len(fips_list)): + pickle_file = os.path.join(pickle_dir, f'{fips_list[0]}_elec_total.pkl') + # Read pickle file + with open(os.path.join(pickle_dir, pickle_file), 'rb') as f: + county_load = pickle.load(f) + county_load = county_load.rename(columns={'total': fips_list[ii]}) + res_elec_load = pd.concat([res_elec_load, county_load], axis=1) + + # Changes in residential load + res_load_change = res_elec_load - res_current_load + res_load_change.index.name = 'Time' + + # Aggregate building load to buses + groupby_dict = dict(zip(fips_list, county_2bus['busIdx'])) + res_load_change_bus = res_load_change.groupby(groupby_dict, axis=1).sum() + res_load_change_bus = res_load_change_bus / 1e3 # convert from kW to MW + + return res_load_change_bus + + def run_nygrid_one_day(s_time, e_time, grid_data, grid_data_dir, opts, init_gen): """ Run NYGrid simulation for one day @@ -105,6 +267,9 @@ def run_nygrid_one_day(s_time, e_time, grid_data, grid_data_dir, opts, init_gen) nygrid_sim.set_gen_ramp_sch(grid_data['genramp30_profile']) nygrid_sim.set_gen_cost_sch(grid_data['gencost0_profile'], grid_data['gencost1_profile']) + if grid_data['genmax_profile_vre'] is not None: + nygrid_sim.set_vre_max_sch(grid_data['genmax_profile_vre']) + # Relax branch flow limits nygrid_sim.relax_external_branch_lim()