diff --git a/data_energy/fitting/data_in_iea_nze.pkl b/data_energy/fitting/data_in_iea_nze.pkl new file mode 100644 index 00000000..3440b57b Binary files /dev/null and b/data_energy/fitting/data_in_iea_nze.pkl differ diff --git a/data_energy/fitting/data_in_iea_nze_after_exec.pkl b/data_energy/fitting/data_in_iea_nze_after_exec.pkl new file mode 100644 index 00000000..3440b57b Binary files /dev/null and b/data_energy/fitting/data_in_iea_nze_after_exec.pkl differ diff --git a/data_energy/fitting/dm_iea_nze.pkl b/data_energy/fitting/dm_iea_nze.pkl new file mode 100644 index 00000000..a25552da Binary files /dev/null and b/data_energy/fitting/dm_iea_nze.pkl differ diff --git a/data_energy/fitting/gaseous_bioenergy.py b/data_energy/fitting/gaseous_bioenergy.py index 74d1750e..b607e6f2 100644 --- a/data_energy/fitting/gaseous_bioenergy.py +++ b/data_energy/fitting/gaseous_bioenergy.py @@ -14,6 +14,8 @@ limitations under the License. ''' import os +import pickle +from copy import deepcopy import numpy as np import pandas as pd @@ -21,9 +23,10 @@ from scipy.interpolate import interp1d from scipy.optimize import minimize from sostrades_core.execution_engine.execution_engine import ExecutionEngine +from sostrades_core.tools.bspline.bspline import BSpline from sostrades_core.tools.post_processing.charts.two_axes_instanciated_chart import ( - InstanciatedSeries, - TwoAxesInstanciatedChart, + InstanciatedSeries, + TwoAxesInstanciatedChart, ) from energy_models.glossaryenergy import GlossaryEnergy @@ -54,106 +57,111 @@ prod_IEA_interpolated = f(years) # increase discretization in order to smooth production between 2020 and 2030 -years_optim = np.arange(years_IEA[0], years_IEA[-1] + 1, 5) #sorted(list(set(years_IEA + list(np.arange(year_start, max(year_start, 2030) + 1))))) +years_optim = np.linspace(year_start, year_end, 8) #np.arange(years_IEA[0], years_IEA[-1] + 1, 5) #sorted(list(set(years_IEA + list(np.arange(year_start, max(year_start, 2030) + 1))))) invest_year_start = 3.432 #G$ -name = 'Test' -model_name = GlossaryEnergy.AnaerobicDigestion +# chose the name so that it mathes the datamanager of the IEA vs NZE study +name = 'usecase_witness_optim_nze_eval' +model_name = f"WITNESS_MDO.WITNESS_Eval.WITNESS.EnergyMix.biogas.{GlossaryEnergy.AnaerobicDigestion}" ns_dict = {'ns_public': name, 'ns_energy': name, 'ns_energy_study': f'{name}', 'ns_biogas': f'{name}', 'ns_resource': name} mod_path = 'energy_models.models.biogas.anaerobic_digestion.anaerobic_digestion_disc.AnaerobicDigestionDiscipline' -ee = ExecutionEngine(name) -ee.ns_manager.add_ns_def(ns_dict) -builder = ee.factory.get_builder_from_module( - model_name, mod_path) -ee.factory.set_builders_to_coupling_builder(builder) +# recover the input data of the discipline from the iea nze scenario +with open('dm_iea_nze.pkl', 'rb') as f: + dm = pickle.load(f) +f.close() +inputs_dict = deepcopy(dm) +inputs_dict.update({f'{name}.{GlossaryEnergy.CO2TaxesValue}': inputs_dict.pop(f'usecase_witness_optim_nze_eval.WITNESS_MDO.WITNESS_Eval.WITNESS.{GlossaryEnergy.CO2TaxesValue}')}) +inputs_dict.update({f'{name}.{GlossaryEnergy.StreamsCO2EmissionsValue}': inputs_dict.pop(f'usecase_witness_optim_nze_eval.WITNESS_MDO.WITNESS_Eval.WITNESS.EnergyMix.{GlossaryEnergy.StreamsCO2EmissionsValue}')}) +inputs_dict.update({f'{name}.{GlossaryEnergy.StreamPricesValue}': inputs_dict.pop(f'usecase_witness_optim_nze_eval.WITNESS_MDO.WITNESS_Eval.WITNESS.EnergyMix.{GlossaryEnergy.StreamPricesValue}')}) +inputs_dict.update({f'{name}.{GlossaryEnergy.ResourcesPriceValue}': inputs_dict.pop(f'usecase_witness_optim_nze_eval.WITNESS_MDO.WITNESS_Eval.WITNESS.EnergyMix.{GlossaryEnergy.ResourcesPriceValue}')}) +inputs_dict.update({f'{name}.{GlossaryEnergy.TransportCostValue}': inputs_dict.pop(f'usecase_witness_optim_nze_eval.WITNESS_MDO.WITNESS_Eval.WITNESS.EnergyMix.biogas.{GlossaryEnergy.TransportCostValue}')}) + + +def run_model(x: list, inputs_dict: dict = inputs_dict, year_end: int = year_end): + invest_before_year_start = x[0:construction_delay] * invest_year_start + invest_years_optim = x[construction_delay:] * invest_year_start + # interpolate on missing years using bspline as in sostrades-core + list_t = np.linspace(0.0, 1.0, len(years)) + bspline = BSpline(n_poles=len(years_optim)) + bspline.set_ctrl_pts(invest_years_optim) + invests, b_array = bspline.eval_list_t(list_t) + invest_df = pd.DataFrame({GlossaryEnergy.Years: years, + GlossaryCore.InvestValue: list(invests)}) -ee.configure() -ee.display_treeview_nodes() + ee = ExecutionEngine(name) + ee.ns_manager.add_ns_def(ns_dict) + builder = ee.factory.get_builder_from_module( + model_name, mod_path) + ee.factory.set_builders_to_coupling_builder(builder) -def run_model(x: list, year_end: int = year_end): - init_prod = x[0] - invest_before_year_start = x[1:1 + construction_delay] - invest_years_optim = x[1 + construction_delay:] - # interpolate on missing years - f = interp1d(years_optim, invest_years_optim, kind='linear') - invests = f(years) - invest_df = pd.DataFrame({GlossaryEnergy.Years: years, - GlossaryCore.InvestValue: list(invests)}) + ee.configure() + #ee.display_treeview_nodes() - inputs_dict = { + inputs_dict.update({ f'{name}.{GlossaryEnergy.YearStart}': year_start, f'{name}.{GlossaryEnergy.YearEnd}': year_end, f'{name}.{model_name}.{GlossaryEnergy.InvestLevelValue}': invest_df, - f'{name}.{GlossaryEnergy.CO2TaxesValue}': pd.DataFrame( - {GlossaryEnergy.Years: years, GlossaryEnergy.CO2Tax: np.linspace(0., 0., len(years))}), - f'{name}.{GlossaryEnergy.StreamsCO2EmissionsValue}': pd.DataFrame({GlossaryEnergy.Years: years, GlossaryEnergy.electricity: np.zeros_like(years), GlossaryEnergy.WetBiomassResource: np.zeros_like(years)}), - f'{name}.{GlossaryEnergy.StreamPricesValue}': pd.DataFrame({GlossaryEnergy.Years: years, GlossaryEnergy.electricity: np.zeros_like(years), GlossaryEnergy.WetBiomassResource: np.zeros_like(years)}), - f'{name}.{GlossaryEnergy.ResourcesPriceValue}': pd.DataFrame({GlossaryEnergy.Years: years, GlossaryEnergy.electricity: np.zeros_like(years), GlossaryEnergy.WetBiomassResource: np.zeros_like(years)}), - f'{name}.{GlossaryEnergy.TransportCostValue}': pd.DataFrame({GlossaryEnergy.Years: years, 'transport': np.zeros(len(years))}), #f'{name}.{model_name}.{GlossaryEnergy.InitialPlantsAgeDistribFactor}': init_age_distrib_factor, - f'{name}.{model_name}.initial_production': init_prod, + f'{name}.{model_name}.initial_production': initial_production, f'{name}.{model_name}.{GlossaryEnergy.InvestmentBeforeYearStartValue}': pd.DataFrame({GlossaryEnergy.Years: np.arange(year_start - construction_delay, year_start), GlossaryEnergy.InvestValue: invest_before_year_start}), - } + }) - # must load the dict twice, otherwise values are not taken into account - ee.load_study_from_input_dict(inputs_dict) ee.load_study_from_input_dict(inputs_dict) ee.execute() prod_df = ee.dm.get_value(ee.dm.get_all_namespaces_from_var_name(GlossaryEnergy.TechnoProductionValue)[0]) #PWh - return prod_df[[GlossaryEnergy.Years, "biogas (TWh)"]], invest_df + return prod_df[[GlossaryEnergy.Years, "biogas (TWh)"]], invest_df, ee def fitting_renewable(x: list): - prod_df, invest_df = run_model(x) + prod_df, invest_df, ee = run_model(x) prod_values_model = prod_df.loc[prod_df[GlossaryEnergy.Years].isin( years_IEA_interpolated), "biogas (TWh)"].values * 1000. # TWh - return (((prod_values_model - prod_IEA_interpolated)) ** 2).mean() + return (((prod_values_model - prod_IEA_interpolated) / (invest_year_start * np.ones_like(prod_values_model))) ** 2).mean() # Initial guess for the variables invest from year 2025 to 2100. -x0 = np.concatenate((np.array([initial_production]), invest_year_start * np.ones(construction_delay), invest_year_start * np.ones(len(years_optim)))) -bounds = [(initial_production * 0.87, initial_production * 0.87)] + [(invest_year_start/2.4, invest_year_start/2.4)] * construction_delay + (len(years_optim)) * [(invest_year_start/3., 3. * invest_year_start)] - +x0 = np.concatenate((np.array([0.]), 1/2.4 * np.ones(construction_delay - 1), np.ones(len(years_optim)))) +bounds = [(0., 0.)] + [(1./2.4/3., 1./2.4 * 3.)] * (construction_delay - 1) + (len(years_optim)) * [(1./3., 3. * 1.)] # Use minimize to find the minimum of the function -result = minimize(fitting_renewable, x0, bounds=bounds, options={'disp': True, 'maxiter': 500, 'maxfun': 500, 'method': 'trust-constr', 'FACTR': 1.e-7}) +result = minimize(fitting_renewable, x0, bounds=bounds, #method='trust-constr', + options={'disp': True, 'maxiter': 500}) #, 'maxfun': 500, 'ftol': 1.e-6, 'maxls': 50}) -prod_df, invest_df = run_model(result.x) +prod_df, invest_df, ee = run_model(result.x) # Print the result print("Function value at the optimum:", result.fun) -print("initial production", result.x[0]) -print("invest before year start", result.x[1:1+construction_delay]) -print("invest at the optimum", result.x[1+construction_delay:]) +print("invest before year start", result.x[0:construction_delay] * invest_year_start) +print("invest at the poles at the optimum", result.x[construction_delay:] * invest_year_start) new_chart = TwoAxesInstanciatedChart('years', 'biogas production (TWh)', - chart_name='Production : model vs historic') + chart_name='Production : witness vs IEA') serie = InstanciatedSeries(list(prod_df[GlossaryEnergy.Years].values), list(prod_df["biogas (TWh)"].values * 1000.), 'model', 'lines+markers') new_chart.series.append(serie) -serie = InstanciatedSeries(years_IEA, df_prod_iea["biogas AnaerobicDigestion (TWh)"].values, 'historic', 'scatter') +serie = InstanciatedSeries(years_IEA, df_prod_iea["biogas AnaerobicDigestion (TWh)"].values, 'IEA', 'scatter') new_chart.series.append(serie) -serie = InstanciatedSeries(list(years_IEA_interpolated), list(prod_IEA_interpolated), 'historic_interpolated', 'lines+markers') +serie = InstanciatedSeries(list(years_IEA_interpolated), list(prod_IEA_interpolated), 'IEA_interpolated', 'lines+markers') new_chart.series.append(serie) new_chart.to_plotly().show() new_chart = TwoAxesInstanciatedChart('years', 'biogas invest (G$)', chart_name='investments') -serie = InstanciatedSeries(list(years_optim), list(result.x)[1+construction_delay:], 'invests_at_poles', 'lines+markers') +serie = InstanciatedSeries(list(years_optim), list(result.x[construction_delay:] * invest_year_start), 'invests_at_poles', 'scatter') new_chart.series.append(serie) -serie = InstanciatedSeries(list(years), list(invest_df[GlossaryEnergy.InvestValue]), 'invests', 'lines') +serie = InstanciatedSeries(list(years), list(invest_df[GlossaryEnergy.InvestValue]), 'invests_bspline', 'lines') new_chart.series.append(serie) new_chart.to_plotly().show() @@ -172,7 +180,5 @@ def fitting_renewable(x: list): df_invest_mix['biogas.AnaerobicDigestion'] = invest_df[GlossaryCore.InvestValue] df_invest_mix.to_csv(invest_mix_csv, index=False, sep=',') # values to set in the invest_design_space_NZE.csv -f = interp1d(years, df_invest_mix['biogas.AnaerobicDigestion'].values, kind='linear') -invest_at_poles = f(np.linspace(year_start, year_end, 8)) -print(f"invest at poles={invest_at_poles}") +print(f"invest at poles={result.x[construction_delay:] * invest_year_start}") diff --git a/data_energy/fitting/hydropower.py b/data_energy/fitting/hydropower.py index eb74e832..9b202a8f 100644 --- a/data_energy/fitting/hydropower.py +++ b/data_energy/fitting/hydropower.py @@ -14,6 +14,8 @@ limitations under the License. ''' import os +import pickle +from copy import deepcopy import numpy as np import pandas as pd @@ -21,9 +23,10 @@ from scipy.interpolate import interp1d from scipy.optimize import minimize from sostrades_core.execution_engine.execution_engine import ExecutionEngine +from sostrades_core.tools.bspline.bspline import BSpline from sostrades_core.tools.post_processing.charts.two_axes_instanciated_chart import ( - InstanciatedSeries, - TwoAxesInstanciatedChart, + InstanciatedSeries, + TwoAxesInstanciatedChart, ) from energy_models.glossaryenergy import GlossaryEnergy @@ -50,111 +53,116 @@ prod_IEA_interpolated = f(years_IEA_interpolated) # increase discretization in order to smooth production between 2020 and 2030 -years_optim = np.arange(years_IEA[0], years_IEA[-1] + 1, 5) #years_IEA_interpolated #sorted(list(set(years_IEA_interpolated + list(np.arange(year_start, max(year_start, 2030) + 1))))) +years_optim = np.linspace(year_start, year_end, 8) #np.arange(years_IEA[0], years_IEA[-1] + 1, 5) #years_IEA_interpolated #sorted(list(set(years_IEA_interpolated + list(np.arange(year_start, max(year_start, 2030) + 1))))) invest_year_start = 18.957 #G$ -name = 'Test' -model_name = GlossaryEnergy.Hydropower -ee = ExecutionEngine(name) -ns_dict = {'ns_public': name, - 'ns_energy': name, - 'ns_energy_study': f'{name}', - 'ns_electricity': name, - 'ns_resource': name} -ee.ns_manager.add_ns_def(ns_dict) +name = 'usecase_witness_optim_nze_eval' +model_name = f"WITNESS_MDO.WITNESS_Eval.WITNESS.EnergyMix.electricity.{GlossaryEnergy.Hydropower}" -mod_path = 'energy_models.models.electricity.hydropower.hydropower_disc.HydropowerDiscipline' -builder = ee.factory.get_builder_from_module( - model_name, mod_path) +# recover the input data of the discipline from the iea nze scenario +with open('dm_iea_nze.pkl', 'rb') as f: + dm = pickle.load(f) +f.close() +inputs_dict = deepcopy(dm) +inputs_dict.update({f'{name}.{GlossaryEnergy.CO2TaxesValue}': inputs_dict.pop(f'usecase_witness_optim_nze_eval.WITNESS_MDO.WITNESS_Eval.WITNESS.{GlossaryEnergy.CO2TaxesValue}')}) +inputs_dict.update({f'{name}.{GlossaryEnergy.StreamsCO2EmissionsValue}': inputs_dict.pop(f'usecase_witness_optim_nze_eval.WITNESS_MDO.WITNESS_Eval.WITNESS.EnergyMix.{GlossaryEnergy.StreamsCO2EmissionsValue}')}) +inputs_dict.update({f'{name}.{GlossaryEnergy.StreamPricesValue}': inputs_dict.pop(f'usecase_witness_optim_nze_eval.WITNESS_MDO.WITNESS_Eval.WITNESS.EnergyMix.{GlossaryEnergy.StreamPricesValue}')}) +inputs_dict.update({f'{name}.{GlossaryEnergy.ResourcesPriceValue}': inputs_dict.pop(f'usecase_witness_optim_nze_eval.WITNESS_MDO.WITNESS_Eval.WITNESS.EnergyMix.{GlossaryEnergy.ResourcesPriceValue}')}) +inputs_dict.update({f'{name}.{GlossaryEnergy.TransportCostValue}': inputs_dict.pop(f'usecase_witness_optim_nze_eval.WITNESS_MDO.WITNESS_Eval.WITNESS.EnergyMix.biogas.{GlossaryEnergy.TransportCostValue}')}) -ee.factory.set_builders_to_coupling_builder(builder) +def run_model(x: list, year_end: int = year_end): + init_prod = x[0] * initial_production + invest_before_year_start = x[1:1 + construction_delay] * invest_year_start + invest_years_optim = x[1 + construction_delay:] * invest_year_start + # interpolate on missing years using bspline as in sostrades-core + list_t = np.linspace(0.0, 1.0, len(years)) + bspline = BSpline(n_poles=len(years_optim)) + bspline.set_ctrl_pts(invest_years_optim) + invests, b_array = bspline.eval_list_t(list_t) + invest_df = pd.DataFrame({GlossaryEnergy.Years: years, + GlossaryCore.InvestValue: list(invests)}) -ee.configure() -ee.display_treeview_nodes() + ee = ExecutionEngine(name) + ns_dict = {'ns_public': name, + 'ns_energy': name, + 'ns_energy_study': f'{name}', + 'ns_electricity': name, + 'ns_resource': name} + ee.ns_manager.add_ns_def(ns_dict) + mod_path = 'energy_models.models.electricity.hydropower.hydropower_disc.HydropowerDiscipline' + builder = ee.factory.get_builder_from_module( + model_name, mod_path) + ee.factory.set_builders_to_coupling_builder(builder) -def run_model(x: list, year_end: int = year_end): - init_prod = x[0] - invest_before_year_start = x[1:1 + construction_delay] - invest_years_optim = x[1 + construction_delay:] - # interpolate on missing years - f = interp1d(years_optim, invest_years_optim, kind='linear') - invests = f(years) - invest_df = pd.DataFrame({GlossaryEnergy.Years: years, - GlossaryCore.InvestValue: list(invests)}) + ee.configure() + #ee.display_treeview_nodes() - inputs_dict = { + inputs_dict.update({ f'{name}.{GlossaryEnergy.YearStart}': year_start, f'{name}.{GlossaryEnergy.YearEnd}': year_end, f'{name}.{model_name}.{GlossaryEnergy.InvestLevelValue}': invest_df, - f'{name}.{GlossaryEnergy.CO2TaxesValue}': pd.DataFrame( - {GlossaryEnergy.Years: years, GlossaryEnergy.CO2Tax: np.linspace(0., 0., len(years))}), - f'{name}.{GlossaryEnergy.StreamsCO2EmissionsValue}': pd.DataFrame({GlossaryEnergy.Years: years}), - f'{name}.{GlossaryEnergy.StreamPricesValue}': pd.DataFrame({GlossaryEnergy.Years: years}), - f'{name}.{GlossaryEnergy.ResourcesPriceValue}': pd.DataFrame({GlossaryEnergy.Years: years}), - f'{name}.{GlossaryEnergy.TransportCostValue}': pd.DataFrame({GlossaryEnergy.Years: years, 'transport': np.zeros(len(years))}), #f'{name}.{model_name}.{GlossaryEnergy.InitialPlantsAgeDistribFactor}': init_age_distrib_factor, f'{name}.{model_name}.initial_production': init_prod, f'{name}.{model_name}.{GlossaryEnergy.InvestmentBeforeYearStartValue}': pd.DataFrame( {GlossaryEnergy.Years: np.arange(year_start - construction_delay, year_start), GlossaryEnergy.InvestValue: invest_before_year_start}), - } - # bug: must load the study twice so that modifications are taked into accout - ee.load_study_from_input_dict(inputs_dict) + }) ee.load_study_from_input_dict(inputs_dict) ee.execute() prod_df = ee.dm.get_value(ee.dm.get_all_namespaces_from_var_name(GlossaryEnergy.TechnoProductionValue)[0]) #PWh - return prod_df[[GlossaryEnergy.Years, "electricity (TWh)"]], invest_df + return prod_df[[GlossaryEnergy.Years, "electricity (TWh)"]], invest_df, ee def fitting_renewable(x: list): - prod_df, invest_df = run_model(x) + prod_df, invest_df, ee = run_model(x) prod_values_model = prod_df.loc[prod_df[GlossaryEnergy.Years].isin( years_IEA_interpolated), "electricity (TWh)"].values * 1000. # TWh - return (((prod_values_model - prod_IEA_interpolated)) ** 2).mean() + return (((prod_values_model - prod_IEA_interpolated) / (initial_production * np.ones_like(prod_values_model))) ** 2).mean() # Initial guess for the variables invest from year 2025 to 2100. -# Initial guess for the variables invest from year 2025 to 2100. -x0 = np.concatenate((np.array([initial_production]), invest_year_start * np.ones(construction_delay), invest_year_start * np.ones(len(years_optim)))) -bounds = [(initial_production, initial_production)] + [(invest_year_start/1., invest_year_start/1.)] * construction_delay + (len(years_optim)) * [(invest_year_start/10., 10. * invest_year_start)] - +# there is a bug with the invest before year start => first value must be set to 0 +# otherwise initial production at year start is not as expected +x0 = np.concatenate((np.array([1.]), np.array([0.]), 80./invest_year_start * np.ones(construction_delay - 1), np.ones(len(years_optim)))) +bounds = [(1., 1.)] + [(0., 0.)] + [(80./invest_year_start/2., 80./invest_year_start * 2.)] * (construction_delay - 1) + (len(years_optim)) * [(1./10., 10.)] # Use minimize to find the minimum of the function -result = minimize(fitting_renewable, x0, bounds=bounds, options={'disp': True, 'maxiter': 500, 'maxfun': 500, 'method': 'trust-constr', 'FACTR': 1.e-7}) +result = minimize(fitting_renewable, x0, bounds=bounds, #method='trust-constr', + options={'disp': True, 'maxiter': 2000, 'xtol': 1e-20}) -prod_df, invest_df = run_model(result.x) +prod_df, invest_df, ee = run_model(result.x) # Print the result print("Function value at the optimum:", result.fun) -print("initial production", result.x[0]) -print("invest before year start", result.x[1:1+construction_delay]) -print("invest at the optimum", result.x[1+construction_delay:]) +print("initial production", result.x[0] * initial_production) +print("invest before year start", result.x[1:1+construction_delay] * invest_year_start) +print("invest at the poles at the optimum", result.x[1+construction_delay:] * invest_year_start) new_chart = TwoAxesInstanciatedChart('years', 'hydropower production (TWh)', - chart_name='Production : model vs historic') + chart_name='witness vs IEA') serie = InstanciatedSeries(list(prod_df[GlossaryEnergy.Years].values), list(prod_df["electricity (TWh)"].values * 1000.), 'model', 'lines') new_chart.series.append(serie) -serie = InstanciatedSeries(years_IEA, df_prod_iea['electricity (TWh)'].values, 'historic', 'scatter') +serie = InstanciatedSeries(years_IEA, df_prod_iea['electricity (TWh)'].values, 'IEA', 'scatter') new_chart.series.append(serie) -serie = InstanciatedSeries(list(years_IEA_interpolated), list(prod_IEA_interpolated), 'historic_interpolated', 'lines+markers') +serie = InstanciatedSeries(list(years_IEA_interpolated), list(prod_IEA_interpolated), 'IEA_interpolated', 'lines+markers') new_chart.series.append(serie) new_chart.to_plotly().show() new_chart = TwoAxesInstanciatedChart('years', 'hydropower invest (G$)', chart_name='investments') -serie = InstanciatedSeries(list(years_optim), list(result.x)[1+construction_delay:], 'invests_at_poles', 'lines+markers') +serie = InstanciatedSeries(list(years_optim), list(result.x[1+construction_delay:] * invest_year_start), 'invests_at_poles', 'scatter') new_chart.series.append(serie) -serie = InstanciatedSeries(list(years), list(invest_df[GlossaryEnergy.InvestValue]), 'invests', 'lines') +serie = InstanciatedSeries(list(years), list(invest_df[GlossaryEnergy.InvestValue]), 'invests_bspline', 'lines') new_chart.series.append(serie) new_chart.to_plotly().show() @@ -176,6 +184,4 @@ def fitting_renewable(x: list): df_invest_mix.to_csv(invest_mix_csv, index=False, sep=',') # values to set in the invest_design_space_NZE.csv -f = interp1d(years, df_invest_mix['electricity.Hydropower'].values, kind='linear') -invest_at_poles = f(np.linspace(year_start, year_end, 8)) -print(f"invest at poles={invest_at_poles}") \ No newline at end of file +print(f"invest at poles={result.x[1+construction_delay:] * invest_year_start}") \ No newline at end of file diff --git a/data_energy/fitting/windpower.py b/data_energy/fitting/windpower.py index 2d3e59aa..102cc573 100644 --- a/data_energy/fitting/windpower.py +++ b/data_energy/fitting/windpower.py @@ -14,6 +14,8 @@ limitations under the License. ''' import os +import pickle +from copy import deepcopy from functools import reduce import numpy as np @@ -22,12 +24,16 @@ from scipy.interpolate import interp1d from scipy.optimize import minimize from sostrades_core.execution_engine.execution_engine import ExecutionEngine +from sostrades_core.tools.bspline.bspline import BSpline from sostrades_core.tools.post_processing.charts.two_axes_instanciated_chart import ( InstanciatedSeries, TwoAxesInstanciatedChart, ) from energy_models.glossaryenergy import GlossaryEnergy +from energy_models.models.electricity.wind_onshore.wind_onshore_disc import ( + WindOnshoreDiscipline, +) """ This script is used to calibrate the windpower invest so that the electricity production matches the IEA NZE scenario @@ -39,6 +45,7 @@ years_IEA = [2020, 2025, 2030, 2035, 2040, 2045, 2050, 2100] years = np.arange(year_start, year_end + 1) + # source: IEA report NZE2021Ch02 models_path_abs = os.path.dirname(os.path.abspath(__file__)).split(os.sep + "models")[0] df_prod_iea = pd.read_csv( @@ -46,19 +53,25 @@ new_row = pd.DataFrame({'years': [2100], 'electricity (TWh)': [35000.]}) df_prod_iea = pd.concat([df_prod_iea, new_row], ignore_index=True) +df_price_iea = pd.read_csv( + os.path.join(models_path_abs, 'models', 'witness-core', 'climateeconomics', 'data', 'IEA_NZE_electricity_Technologies_Mix_prices.csv')) + # interpolate data between 2050 and 2100 years_IEA_interpolated = years f = interp1d(years_IEA, df_prod_iea['electricity (TWh)'].values, kind='linear') prod_IEA_interpolated = f(years_IEA_interpolated) -# increase discretization in order to smooth production between 2020 and 2030 -years_optim = years_IEA_interpolated #sorted(list(set(years_IEA + list(np.arange(year_start, max(year_start, 2030) + 1))))) +# optimization at the poles just like in witness-full study +years_optim = np.linspace(year_start, year_end, 8) #years_IEA_interpolated #sorted(list(set(years_IEA + list(np.arange(year_start, max(year_start, 2030) + 1))))) invest_year_start = 80. #G$ +construction_delay = GlossaryEnergy.TechnoConstructionDelayDict['WindOffshore'] # same construction delay for windonshore and windoffshore +if construction_delay != GlossaryEnergy.TechnoConstructionDelayDict['WindOnshore']: + raise ValueError("must adapt script as construction delay for windOnshore and windOffshore differ") -name = 'Test' -model_name_onshore = GlossaryEnergy.WindOnshore -model_name_offshore = GlossaryEnergy.WindOffshore +name = 'usecase_witness_optim_nze_eval' +model_name_onshore = f"WITNESS_MDO.WITNESS_Eval.WITNESS.EnergyMix.electricity.{GlossaryEnergy.WindOnshore}" +model_name_offshore = f"WITNESS_MDO.WITNESS_Eval.WITNESS.EnergyMix.electricity.{GlossaryEnergy.WindOffshore}" ns_dict = {'ns_public': name, 'ns_energy': name, 'ns_energy_study': f'{name}', @@ -67,57 +80,76 @@ mod_path_onshore = 'energy_models.models.electricity.wind_onshore.wind_onshore_disc.WindOnshoreDiscipline' mod_path_offshore = 'energy_models.models.electricity.wind_offshore.wind_offshore_disc.WindOffshoreDiscipline' -ee = ExecutionEngine(name) -ee.ns_manager.add_ns_def(ns_dict) -builder = [] -builder.append(ee.factory.get_builder_from_module( - model_name_onshore, mod_path_onshore)) -builder.append(ee.factory.get_builder_from_module( - model_name_offshore, mod_path_offshore)) -ee.factory.set_builders_to_coupling_builder(builder) - -ee.configure() -ee.display_treeview_nodes() - - - +# if want to modify the capex of both onshore and offshore +#dict_techno_dict_default = {model_name_onshore: WindOnshoreDiscipline.techno_infos_dict_default, +# model_name_offshore: WindOffshoreDiscipline.techno_infos_dict_default} +techno_info_dict_default = WindOnshoreDiscipline.techno_infos_dict_default +Capex_init0 = WindOnshoreDiscipline.techno_infos_dict_default['Capex_init'] + +# recover the input data of the discipline from the iea nze scenario +with open('dm_iea_nze.pkl', 'rb') as f: + dm = pickle.load(f) +f.close() +inputs_dict = deepcopy(dm) +inputs_dict.update({f'{name}.{GlossaryEnergy.CO2TaxesValue}': inputs_dict.pop(f'usecase_witness_optim_nze_eval.WITNESS_MDO.WITNESS_Eval.WITNESS.{GlossaryEnergy.CO2TaxesValue}')}) +inputs_dict.update({f'{name}.{GlossaryEnergy.StreamsCO2EmissionsValue}': inputs_dict.pop(f'usecase_witness_optim_nze_eval.WITNESS_MDO.WITNESS_Eval.WITNESS.EnergyMix.{GlossaryEnergy.StreamsCO2EmissionsValue}')}) +inputs_dict.update({f'{name}.{GlossaryEnergy.StreamPricesValue}': inputs_dict.pop(f'usecase_witness_optim_nze_eval.WITNESS_MDO.WITNESS_Eval.WITNESS.EnergyMix.{GlossaryEnergy.StreamPricesValue}')}) +inputs_dict.update({f'{name}.{GlossaryEnergy.ResourcesPriceValue}': inputs_dict.pop(f'usecase_witness_optim_nze_eval.WITNESS_MDO.WITNESS_Eval.WITNESS.EnergyMix.{GlossaryEnergy.ResourcesPriceValue}')}) +inputs_dict.update({f'{name}.{GlossaryEnergy.TransportCostValue}': inputs_dict.pop(f'usecase_witness_optim_nze_eval.WITNESS_MDO.WITNESS_Eval.WITNESS.EnergyMix.electricity.{GlossaryEnergy.TransportCostValue}')}) + +# initial_production from IEA is split between windonshore and windoffshore following arbitrary ratio +init_prod_onshore_over_offshore = 1508. / 107.69 # taken from initial witness results +initial_prod = df_prod_iea.loc[df_prod_iea[GlossaryEnergy.Years] == year_start]['electricity (TWh)'].values[0] +init_prod_dict = { + model_name_onshore: initial_prod * init_prod_onshore_over_offshore / (1. + init_prod_onshore_over_offshore), + model_name_offshore: initial_prod / (1. + init_prod_onshore_over_offshore)} +ratio_invest_onshore_offshore = 3.6689 # taken from initial witness results def run_model(x: list, year_end: int = year_end): - init_age_distrib_factor = x[0] - invest_years_optim = x[1:] - # interpolate on missing years - f = interp1d(years_optim, invest_years_optim, kind='linear') - invests = f(years) + techno_info_dict_default['Capex_init'] = Capex_init0 * x[0] + invest_before_year_start = x[1:construction_delay + 1] * invest_year_start + invest_years_optim = x[construction_delay + 1:] * invest_year_start + # interpolate on missing years using bspline as in sostrades-core + list_t = np.linspace(0.0, 1.0, len(years)) + bspline = BSpline(n_poles=len(years_optim)) + bspline.set_ctrl_pts(invest_years_optim) + invests, b_array = bspline.eval_list_t(list_t) invest_df = pd.DataFrame({GlossaryEnergy.Years: years, GlossaryCore.InvestValue: list(invests)}) - + invest_before_year_start_df = pd.DataFrame({GlossaryEnergy.Years: np.arange(year_start - construction_delay, year_start), + GlossaryEnergy.InvestValue: invest_before_year_start}) # split investment between onshore and offshore - ratio_invest_onshore_offshore = 3.6689 * np.ones(len(invests)) # taken from initial witness results - invest_df[model_name_onshore] = invest_df[GlossaryCore.InvestValue] * ratio_invest_onshore_offshore / (1. + ratio_invest_onshore_offshore) - invest_df[model_name_offshore] = invest_df[GlossaryCore.InvestValue] / (1. + ratio_invest_onshore_offshore) - init_prod_onshore_over_offshore = 1508./107.69 # taken from initial witness results - initial_prod = df_prod_iea.loc[df_prod_iea[GlossaryEnergy.Years] == year_start]['electricity (TWh)'].values[0] - init_prod_dict = {model_name_onshore: initial_prod * init_prod_onshore_over_offshore / (1. + init_prod_onshore_over_offshore), - model_name_offshore: initial_prod / (1. + init_prod_onshore_over_offshore)} - - inputs_dict = { + for df in [invest_df, invest_before_year_start_df]: + df[model_name_onshore] = df[GlossaryCore.InvestValue] * ratio_invest_onshore_offshore / (1. + ratio_invest_onshore_offshore) + df[model_name_offshore] = df[GlossaryCore.InvestValue] / (1. + ratio_invest_onshore_offshore) + + + ee = ExecutionEngine(name) + ee.ns_manager.add_ns_def(ns_dict) + builder = [] + builder.append(ee.factory.get_builder_from_module( + model_name_onshore, mod_path_onshore)) + builder.append(ee.factory.get_builder_from_module( + model_name_offshore, mod_path_offshore)) + ee.factory.set_builders_to_coupling_builder(builder) + + ee.configure() + #ee.display_treeview_nodes() + + inputs_dict.update({ f'{name}.{GlossaryEnergy.YearStart}': year_start, f'{name}.{GlossaryEnergy.YearEnd}': year_end, - f'{name}.{GlossaryEnergy.CO2TaxesValue}': pd.DataFrame( - {GlossaryEnergy.Years: years, GlossaryEnergy.CO2Tax: np.linspace(0., 0., len(years))}), - f'{name}.{GlossaryEnergy.StreamsCO2EmissionsValue}': pd.DataFrame({GlossaryEnergy.Years: years}), - f'{name}.{GlossaryEnergy.StreamPricesValue}': pd.DataFrame({GlossaryEnergy.Years: years}), - f'{name}.{GlossaryEnergy.ResourcesPriceValue}': pd.DataFrame({GlossaryEnergy.Years: years}), - f'{name}.{GlossaryEnergy.TransportCostValue}': pd.DataFrame({GlossaryEnergy.Years: years, 'transport': np.zeros(len(years))}), - } + f'{name}.{model_name_onshore}.techno_infos_dict': techno_info_dict_default, + }) for model_name in [model_name_offshore, model_name_onshore]: inputs_dict.update({ - f'{name}.{model_name}.{GlossaryEnergy.InitialPlantsAgeDistribFactor}': init_age_distrib_factor, + #f'{name}.{model_name}.{GlossaryEnergy.InitialPlantsAgeDistribFactor}': init_age_distrib_factor, f'{name}.{model_name}.initial_production': init_prod_dict[model_name], f'{name}.{model_name}.{GlossaryEnergy.InvestLevelValue}': pd.DataFrame({GlossaryEnergy.Years: years, GlossaryCore.InvestValue: invest_df[model_name].values}), + f'{name}.{model_name}.{GlossaryEnergy.InvestmentBeforeYearStartValue}': pd.DataFrame({GlossaryEnergy.Years: np.arange(year_start - construction_delay, year_start), + GlossaryEnergy.InvestValue:invest_before_year_start_df[model_name].values}), }) - # bug: must load the study twice so that modifications are taked into accout - ee.load_study_from_input_dict(inputs_dict) + ee.load_study_from_input_dict(inputs_dict) ee.execute() @@ -130,49 +162,66 @@ def run_model(x: list, year_end: int = year_end): df_prod['electricity (TWh)'] = df_prod.drop(GlossaryEnergy.Years, axis=1).sum(axis=1) * 1000. #PWh df_prod_model = df_prod.loc[df_prod[GlossaryEnergy.Years].isin(years_IEA_interpolated)] - return df_prod, df_prod_model, invest_df + price_df = ee.dm.get_value(f"{name}.{model_name}.{GlossaryEnergy.TechnoPricesValue}") + + return df_prod, price_df, df_prod_model, invest_df, ee def fitting_renewable(x: list): - df_prod, df_prod_model, invest_df = run_model(x) - return (((df_prod_model['electricity (TWh)'].values - prod_IEA_interpolated)) ** 2).mean() + df_prod, price_df, df_prod_model, invest_df, ee = run_model(x) + price_iea_values = df_price_iea['WindOnshore'].values + years_price_iea = df_price_iea['years'].values + price_model_values = (price_df.loc[price_df[GlossaryEnergy.Years].isin(years_price_iea), f"{GlossaryEnergy.WindOnshore}_wotaxes"]).values + + return ((((df_prod_model['electricity (TWh)'].values - prod_IEA_interpolated)/prod_IEA_interpolated.mean()) ** 2).mean() + (((price_model_values - price_iea_values)/price_iea_values.mean()) ** 2).mean()) # Initial guess for the variables invest from year 2025 to 2100. -x0 = np.concatenate((np.array([1.0]), invest_year_start * np.ones(len(years_optim)))) -bounds = [(1., 2.)] + (len(years_optim)) * [(invest_year_start/10., 10.0 * invest_year_start)] +x0 = np.concatenate((np.array([1.]), np.array([0.]), 80./invest_year_start * np.ones(construction_delay - 1), np.ones(len(years_optim)))) +bounds = [(0.5, 1.5)] + [(0., 0.)] + [(80./invest_year_start/2., 80./invest_year_start * 2.)] * (construction_delay - 1) + (len(years_optim)) * [(1./10., 10.)] # Use minimize to find the minimum of the function result = minimize(fitting_renewable, x0, bounds=bounds) -df_prod, df_prod_model, invest_df = run_model(result.x) +df_prod, price_df, df_prod_model, invest_df, ee = run_model(result.x) # Print the result print("Function value at the optimum:", result.fun) -print("init age distrib at the optimum", result.x[0]) -print("invest at the optimum", result.x[1:]) -print("prod at the optimum", df_prod_model['electricity (TWh)'].values) +print('Capex_init wind onshore:', result.x[0] * Capex_init0) +print("invest before year start", result.x[1:construction_delay + 1] * invest_year_start) +print("invest at the poles at the optimum", result.x[construction_delay + 1:] * invest_year_start) + new_chart = TwoAxesInstanciatedChart('years', 'production (TWh)', - chart_name='Windpower Production : model vs historic') + chart_name='Windpower Production : witness vs IEA') -serie = InstanciatedSeries(years_IEA_interpolated, df_prod_model['electricity (TWh)'].values, 'model', 'lines') +serie = InstanciatedSeries(list(years_IEA_interpolated), df_prod_model['electricity (TWh)'].values, 'model', 'lines') new_chart.series.append(serie) -serie = InstanciatedSeries(years_IEA, df_prod_iea['electricity (TWh)'].values, 'historic', 'scatter') +serie = InstanciatedSeries(years_IEA, df_prod_iea['electricity (TWh)'].values, 'IEA', 'scatter') new_chart.series.append(serie) -serie = InstanciatedSeries(list(years_IEA_interpolated), list(prod_IEA_interpolated), 'historic_interpolated', 'lines+markers') +serie = InstanciatedSeries(list(years_IEA_interpolated), list(prod_IEA_interpolated), 'IEA_interpolated', 'lines+markers') new_chart.series.append(serie) new_chart.to_plotly().show() new_chart = TwoAxesInstanciatedChart('years', 'invest (G$)', chart_name='Windpower investments') -serie = InstanciatedSeries(years_optim, list(result.x)[1:], 'invests_at_poles', 'lines+markers') +serie = InstanciatedSeries(list(years_optim), list(result.x[construction_delay + 1:] * invest_year_start), 'invests_at_poles', 'scatter') +new_chart.series.append(serie) +serie = InstanciatedSeries(list(years), list(invest_df[GlossaryEnergy.InvestValue]), 'invests', 'lines') +new_chart.series.append(serie) + +new_chart.to_plotly().show() + +new_chart = TwoAxesInstanciatedChart('years', 'Price ($/MWh)', + chart_name='Wind Onshore price') +serie = InstanciatedSeries(list(df_price_iea['years'].values), list(df_price_iea['WindOnshore'].values), 'IEA', 'scatter') new_chart.series.append(serie) -serie = InstanciatedSeries(years, list(invest_df[GlossaryEnergy.InvestValue]), 'invests', 'lines') +# in witness vs iea post -processing, take f"{GlossaryEnergy.WindOnshore}" but same value +serie = InstanciatedSeries(list(years), list(price_df[f"{GlossaryEnergy.WindOnshore}_wotaxes"].values), 'Witness', 'lines') new_chart.series.append(serie) new_chart.to_plotly().show() @@ -198,7 +247,7 @@ def fitting_renewable(x: list): df_invest_mix['electricity.WindOnshore'] = invest_df['WindOnshore'] df_invest_mix.to_csv(invest_mix_csv, index=False, sep=',') # values to set in the invest_design_space_NZE.csv -for techno in ['WindOffshore', 'WindOnshore']: - f = interp1d(years, df_invest_mix[f"electricity.{techno}"].values, kind='linear') - invest_at_poles = f(np.linspace(year_start, year_end, 8)) - print(f"invest at poles for {techno}={invest_at_poles}") \ No newline at end of file +print(f"invest at poles for WindOnshore={result.x[construction_delay + 1:] * invest_year_start * ratio_invest_onshore_offshore / (1. + ratio_invest_onshore_offshore)}") +print(f"invest at poles for WindOffshore={result.x[construction_delay + 1:] * invest_year_start / (1. + ratio_invest_onshore_offshore)}") +print(f"invest before year start for WindOnshore={result.x[1:construction_delay +1] * invest_year_start * ratio_invest_onshore_offshore / (1. + ratio_invest_onshore_offshore)}") +print(f"invest before year start for WindOffshore={result.x[1:construction_delay + 1] * invest_year_start / (1. + ratio_invest_onshore_offshore)}") \ No newline at end of file diff --git a/energy_models/models/liquid_fuel/refinery/refinery.py b/energy_models/models/liquid_fuel/refinery/refinery.py index e573b7c3..0a871976 100644 --- a/energy_models/models/liquid_fuel/refinery/refinery.py +++ b/energy_models/models/liquid_fuel/refinery/refinery.py @@ -40,6 +40,24 @@ def __init__(self, name): self.dprod_dinvest = None self.dprod_list_dcapex_list = None + def compute_cost_of_resources_usage(self): + """ + Cost of resource R = need of resource R x price of resource R + + Does not take oil price into account + """ + cost_of_resource_usage = { + GlossaryEnergy.Years: self.years, + } + for resource in self.resources_used_for_production: + if resource == GlossaryEnergy.OilResource: + # Skip OilResource so not to count it twice + cost_of_resource_usage[resource] = 0.0 + else: + cost_of_resource_usage[resource] = self.cost_details[f"{resource}_needs"].values * self.resources_prices[resource].values + + self.cost_of_resources_usage = pd.DataFrame(cost_of_resource_usage) + def get_fuel_needs(self): """ Get the fuel needs for 1 kwh of the energy producted by the technology diff --git a/energy_models/models/methane/fossil_gas/fossil_gas.py b/energy_models/models/methane/fossil_gas/fossil_gas.py index 03fb8ecd..e395c384 100644 --- a/energy_models/models/methane/fossil_gas/fossil_gas.py +++ b/energy_models/models/methane/fossil_gas/fossil_gas.py @@ -15,6 +15,8 @@ limitations under the License. ''' +import pandas as pd + from energy_models.core.stream_type.carbon_models.carbon_capture import CarbonCapture from energy_models.core.stream_type.energy_models.methane import Methane from energy_models.core.techno_type.base_techno_models.methane_techno import ( @@ -39,6 +41,24 @@ def get_fuel_needs(self): return fuel_need + def compute_cost_of_resources_usage(self): + """ + Cost of resource R = need of resource R x price of resource R + + Does not take natural gas price into account + """ + cost_of_resource_usage = { + GlossaryEnergy.Years: self.years, + } + for resource in self.resources_used_for_production: + if resource == GlossaryEnergy.NaturalGasResource: + # Skip NaturalGasResource so not to count it twice + cost_of_resource_usage[resource] = 0.0 + else: + cost_of_resource_usage[resource] = self.cost_details[f"{resource}_needs"].values * self.resources_prices[resource].values + + self.cost_of_resources_usage = pd.DataFrame(cost_of_resource_usage) + def compute_resources_needs(self): self.cost_details[f'{self.NATURAL_GAS_RESOURCE_NAME}_needs'] = self.get_fuel_needs() / Methane.data_energy_dict['calorific_value'] # kg/kWh diff --git a/platform_version_required.txt b/platform_version_required.txt index 8a58a0dc..15a2b33b 100644 --- a/platform_version_required.txt +++ b/platform_version_required.txt @@ -1 +1 @@ -v4.1.3 \ No newline at end of file +v4.2.0 \ No newline at end of file