diff --git a/MINDSET_FTT_Power/SourceCode/Power/ftt_p_main.py b/MINDSET_FTT_Power/SourceCode/Power/ftt_p_main.py index 09a0f42..0e17346 100644 --- a/MINDSET_FTT_Power/SourceCode/Power/ftt_p_main.py +++ b/MINDSET_FTT_Power/SourceCode/Power/ftt_p_main.py @@ -56,6 +56,10 @@ Main solution function for the module """ +# Standard library imports +import configparser +from pathlib import Path + # Third party imports import numpy as np @@ -73,6 +77,25 @@ +# Settings parsed once at import time ÔÇö not per solve() call +_PACKAGE_ROOT = Path(__file__).parents[2] +_config = configparser.ConfigParser() +_config.read(str(_PACKAGE_ROOT / "settings.ini")) + +_PRSC_BASE_YEAR = int(_config.get('settings', 'prsc_base_year', fallback='2013')) +_EX_BASE_YEAR = int(_config.get('settings', 'ex_base_year', fallback='2018')) +_PRSC_DEFL_YEAR = int(_config.get('settings', 'prsc_defl_year', fallback='2015')) +_RLDC_START_YEAR = int(_config.get('settings', 'rldc_start_year', fallback='2018')) +_BCET_COPY_RANGE_END = int(_config.get('settings', 'bcet_copy_range_end', fallback='17')) +_USD_EUR_RATE = float(_config.get('settings', 'usd_eur_rate', fallback='1.0018')) + +_PRSC_VAR = f"PRSC{str(_PRSC_BASE_YEAR)[2:]}" # e.g. "PRSC13" +_PRSC_EX_VAR = f"PRSC{str(_EX_BASE_YEAR)[2:]}" # e.g. "PRSC18" +_EX_VAR = f"EX{str(_EX_BASE_YEAR)[2:]}" # e.g. "EX18" +_REX_VAR = f"REX{str(_EX_BASE_YEAR)[2:]}" # e.g. "REX18" +_PRSC_DEFL_VAR = f"PRSC{str(_PRSC_DEFL_YEAR)[2:]}" # e.g. "PRSC15" + + # %% main function # ----------------------------------------------------------------------------- # ----------------------------- Main ------------------------------------------ @@ -141,10 +164,10 @@ def solve(data, time_lag, iter_lag, titles, conv, histend, year, domain): # Copy over PRSC/EX values - data['PRSC18'] = np.copy(time_lag['PRSC18'] ) - data['EX18'] = np.copy(time_lag['EX18'] ) - data['PRSC15'] = np.copy(time_lag['PRSC15'] ) - data["REX18"] = np.copy(time_lag["REX18"]) + data[_PRSC_EX_VAR] = np.copy(time_lag[_PRSC_EX_VAR]) + data[_EX_VAR] = np.copy(time_lag[_EX_VAR]) + data[_PRSC_DEFL_VAR] = np.copy(time_lag[_PRSC_DEFL_VAR]) + data[_REX_VAR] = np.copy(time_lag[_REX_VAR]) data['FPIX'][data['FPIX'] == 0] = 1 time_lag['FPIX'][time_lag['FPIX'] == 0] = 1 # Increase fuel prices @@ -155,13 +178,13 @@ def solve(data, time_lag, iter_lag, titles, conv, histend, year, domain): T_Scal = 10 # Time scaling factor used in the share dynamics # Initialisation, which corresponds to lines 389 to 556 in fortran - if year == 2013: - data['PRSC13'] = np.copy(data['PRSCX']) + if year == _PRSC_BASE_YEAR: + data[_PRSC_VAR] = np.copy(data['PRSCX']) - if year == 2018: - data['PRSC18'] = np.copy(data['PRSCX']) - data['EX18'] = np.copy(data['EXX']) - data['REX18'] = np.copy(data['REXX']) + if year == _EX_BASE_YEAR: + data[_PRSC_EX_VAR] = np.copy(data['PRSCX']) + data[_EX_VAR] = np.copy(data['EXX']) + data[_REX_VAR] = np.copy(data['REXX']) data['MEWL'][:, :, 0] = data["MWLO"][:, :, 0] data['MEWK'][:, :, 0] = np.divide(data['MEWG'][:, :, 0], data['MEWL'][:, :, 0], @@ -255,15 +278,12 @@ def solve(data, time_lag, iter_lag, titles, conv, histend, year, domain): #%% # Up to the last year of historical market share data elif year <= histend['MEWG']: - if year == 2015: - data['PRSC15'] = np.copy(data['PRSCX']) + if year == _PRSC_DEFL_YEAR: + data[_PRSC_DEFL_VAR] = np.copy(data['PRSCX']) - # Set starting values for MERC - data['MERC'][:, 0, 0] = 0.255 - data['MERC'][:, 1, 0] = 5.689 - data['MERC'][:, 2, 0] = 0.4246 - data['MERC'][:, 3, 0] = 3.374 + # Set starting values for MERC from lagged data (avoids hardcoded floats) + data['MERC'][:, :4, 0] = time_lag['MERC'][:, :4, 0] data['MERC'][:, 4, 0] = 0.001 data['MERC'][:, 7, 0] = 0.001 # Calculate electricty trade shares @@ -277,7 +297,7 @@ def solve(data, time_lag, iter_lag, titles, conv, histend, year, domain): # loadfac = data['MWLO'][:, :, 0] # data['MEWL'][:, :, 0] = np.copy(loadfac) - if year > 2018: + if year > _EX_BASE_YEAR: data['MEWL'][:, :, 0] = time_lag['MEWL'][:, :, 0].copy() cond = np.logical_and(data['MEWL'][:, :, 0] < 0.01, data['MWLO'][:, :, 0] > 0.0) @@ -300,14 +320,14 @@ def solve(data, time_lag, iter_lag, titles, conv, histend, year, domain): # Call RLDC function for capacity and load factor by LB, and storage costs - if year >= 2018: + if year >= _RLDC_START_YEAR: # 1 and 2 -- Estimate RLDC and storage parameters data = rldc(data, time_lag, iter_lag, year, titles) # 3--- Call dispatch routine to connect market shares to load bands # Call DSPCH function to dispatch flexible capacity based on MC - if year == 2018: + if year == _RLDC_START_YEAR: mslb, mllb, mes1, mes2 = dspch(data['MWDD'], data['MEWS'], data['MKLB'], data['MCRT'], data['MEWL'], data['MWMC'], data['MMCD'], len(titles['RTI']), len(titles['T2TI']), len(titles['LBTI'])) @@ -321,12 +341,12 @@ def solve(data, time_lag, iter_lag, titles, conv, histend, year, domain): data['MES2'] = mes2 # Change currency from EUR2015 to USD2013 - if year >= 2015: - # usa_idx = titles['RTI_short'].index('USA') - data['MSSP'][:, :, 0] = data['MSSP'][:, :, 0] * (data['PRSC13'][:, 0, 0, np.newaxis]/data['PRSC15'][:, 0, 0, np.newaxis])/ 1.0018 - data['MLSP'][:, :, 0] = data['MLSP'][:, :, 0] * (data['PRSC13'][:, 0, 0, np.newaxis]/data['PRSC15'][:, 0, 0, np.newaxis])/ 1.0018 - data['MSSM'][:, :, 0] = data['MSSM'][:, :, 0] * (data['PRSC13'][:, 0, 0, np.newaxis]/data['PRSC15'][:, 0, 0, np.newaxis])/ 1.0018 - data['MLSM'][:, :, 0] = data['MLSM'][:, :, 0] * (data['PRSC13'][:, 0, 0, np.newaxis]/data['PRSC15'][:, 0, 0, np.newaxis])/ 1.0018 + if year >= _PRSC_DEFL_YEAR: + _prsc_ratio = data[_PRSC_VAR][:, 0, 0, np.newaxis] / data[_PRSC_DEFL_VAR][:, 0, 0, np.newaxis] + data['MSSP'][:, :, 0] = data['MSSP'][:, :, 0] * _prsc_ratio / _USD_EUR_RATE + data['MLSP'][:, :, 0] = data['MLSP'][:, :, 0] * _prsc_ratio / _USD_EUR_RATE + data['MSSM'][:, :, 0] = data['MSSM'][:, :, 0] * _prsc_ratio / _USD_EUR_RATE + data['MLSM'][:, :, 0] = data['MLSM'][:, :, 0] * _prsc_ratio / _USD_EUR_RATE # TODO: This is not per se correct but it's how it is in E3ME else: @@ -457,7 +477,7 @@ def solve(data, time_lag, iter_lag, titles, conv, histend, year, domain): data["MEWW"][0, :, 0] = time_lag['MEWW'][0, :, 0] + dw # Copy over the technology cost categories that do not change (all except prices which are updated through learning-by-doing below) - data['BCET'][:, :, 1:17] = time_lag['BCET'][:, :, 1:17].copy() + data['BCET'][:, :, 1:_BCET_COPY_RANGE_END] = time_lag['BCET'][:, :, 1:_BCET_COPY_RANGE_END].copy() # Store gamma values in the cost matrix (in case it varies over time) data['BCET'][:, :, c2ti['21 Gamma ($/MWh)']] = data['MGAM'][:, :, 0] @@ -641,13 +661,12 @@ def solve(data, time_lag, iter_lag, titles, conv, histend, year, domain): # Call RLDC function for capacity and load factor by LB, and storage costs data = rldc(data, time_lag, data_dt, year, titles) - # Change currency from EUR2015 to USD2013 (This is wrong, but in terms of logic and by misstating currency year for storage) - # usa_idx = titles['RTI_short'].index('USA') - - data['MSSP'][:, :, 0] = data['MSSP'][:, :, 0] * (data['PRSC13'][:, 0, 0, np.newaxis]/data['PRSC15'][:, 0, 0, np.newaxis]) / 1.0018 - data['MLSP'][:, :, 0] = data['MLSP'][:, :, 0] * (data['PRSC13'][:, 0, 0, np.newaxis]/data['PRSC15'][:, 0, 0, np.newaxis]) / 1.0018 - data['MSSM'][:, :, 0] = data['MSSM'][:, :, 0] * (data['PRSC13'][:, 0, 0, np.newaxis]/data['PRSC15'][:, 0, 0, np.newaxis]) / 1.0018 - data['MLSM'][:, :, 0] = data['MLSM'][:, :, 0] * (data['PRSC13'][:, 0, 0, np.newaxis]/data['PRSC15'][:, 0, 0, np.newaxis]) / 1.0018 + # Change currency from EUR2015 to USD2013 (acknowledged approximation) + _prsc_ratio = data[_PRSC_VAR][:, 0, 0, np.newaxis] / data[_PRSC_DEFL_VAR][:, 0, 0, np.newaxis] + data['MSSP'][:, :, 0] = data['MSSP'][:, :, 0] * _prsc_ratio / _USD_EUR_RATE + data['MLSP'][:, :, 0] = data['MLSP'][:, :, 0] * _prsc_ratio / _USD_EUR_RATE + data['MSSM'][:, :, 0] = data['MSSM'][:, :, 0] * _prsc_ratio / _USD_EUR_RATE + data['MLSM'][:, :, 0] = data['MLSM'][:, :, 0] * _prsc_ratio / _USD_EUR_RATE # ================================================================= @@ -788,7 +807,7 @@ def solve(data, time_lag, iter_lag, titles, conv, histend, year, domain): # Copy over the technology cost categories. We update the investment and capacity factors below - data['BCET'][:, :, 1:17] = time_lag['BCET'][:, :, 1:17].copy() + data['BCET'][:, :, 1:_BCET_COPY_RANGE_END] = time_lag['BCET'][:, :, 1:_BCET_COPY_RANGE_END].copy() # Store gamma values in the cost matrix (in case it varies over time) data['BCET'][:, :, c2ti['21 Gamma ($/MWh)']] = data['MGAM'][:, :, 0] diff --git a/MINDSET_FTT_Power/SourceCode/model_class.py b/MINDSET_FTT_Power/SourceCode/model_class.py index 43ea299..6422598 100644 --- a/MINDSET_FTT_Power/SourceCode/model_class.py +++ b/MINDSET_FTT_Power/SourceCode/model_class.py @@ -13,16 +13,14 @@ # Standard library imports import configparser -import copy +from pathlib import Path # Third party imports import numpy as np from tqdm import tqdm -# Local library imports -# Separate FTT modules -import MINDSET_FTT_Power.SourceCode.Power.ftt_p_main as ftt_p - +from ftt_source.paths import set_paths +from ftt_source.Power.ftt_p_main import solve as ftt_p_solve, build_power_settings # Support modules import MINDSET_FTT_Power.SourceCode.support.input_functions as in_f @@ -94,8 +92,10 @@ def __init__(self): """ Instantiate model run object """ # Attributes given in settings.ini file + _mindset_root = Path(__file__).parents[1] + _settings_ini = str(_mindset_root / 'settings.ini') config = configparser.ConfigParser() - config.read('MINDSET_FTT_Power/settings.ini') + config.read(_settings_ini) self.name = config.get('settings', 'name') self.model_start = int(config.get('settings', 'model_start')) self.model_end = int(config.get('settings', 'model_end')) @@ -107,21 +107,47 @@ def __init__(self): self.ftt_modules = config.get('settings', 'enable_modules') self.scenarios = config.get('settings', 'scenarios') + # Point FTT_Standalone at MINDSET's Utilities and settings (once at startup) + set_paths(utilities_path=str(_mindset_root / 'Utilities')) + self.ftt_settings_path = _settings_ini + self.carbon_price_conv_factor = float( + config.get('settings', 'carbon_price_conv_factor', fallback='1.3281')) + # Read exchange-rate variable names from settings (must match ftt_p_main.py) + _ex_base_year = int(config.get('settings', 'ex_base_year', fallback='2018')) + self._prsc_ex_var = f"PRSC{str(_ex_base_year)[2:]}" # e.g. "PRSC18" + self._ex_var = f"EX{str(_ex_base_year)[2:]}" # e.g. "EX18" + # Load classification titles self.titles = titles_f.load_titles() self.conv = titles_f.load_converters() + self.power_settings = build_power_settings(self.titles, config) # Load variable dimensions self.dims, self.histend, self.domain, self.forstart = dims_f.load_dims() - - # # Set up csv files if they do not exist yet - # initialise_csv_files(self.ftt_modules, self.scenarios) - - # Retrieve inputs + + # Retrieve inputs ÔÇö C2TI size must match the CSV files at this point self.input = in_f.load_data(self.titles, self.dims, self.timeline, self.scenarios, self.ftt_modules, self.forstart) + # After loading, extend BCET and C2TI so FTT_Standalone's ftt_p_lcoe can find + # '22 Gamma' (index 21) and '23 Value factor' (index 22) by name. + # The CSV data ends at '21 Gamma ($/MWh)' (index 20); we copy it to column 21 + # and default Value factor to 1.0 for all technologies. + c2ti_list = list(self.titles['C2TI']) + gamma_col = c2ti_list.index('21 Gamma ($/MWh)') # index 20 + new_entries = [e for e in ('22 Gamma', '23 Value factor') if e not in c2ti_list] + if new_entries: + n_extra = len(new_entries) + for scen in self.input: + bcet = self.input[scen]['BCET'] # (RTI, T2TI, C2TI, 1) + extra = np.ones((bcet.shape[0], bcet.shape[1], n_extra, bcet.shape[3])) + extra[:, :, 0, :] = bcet[:, :, gamma_col, :] # '22 Gamma' = copy of col 21 + # '23 Value factor' stays 1.0 + self.input[scen]['BCET'] = np.concatenate([bcet, extra], axis=2) + c2ti_list.extend(new_entries) + self.titles['C2TI'] = tuple(c2ti_list) + # Initialize remaining attributes self.variables = {} @@ -130,6 +156,17 @@ def __init__(self): self.output = {scen: {var: np.full_like(self.input[scen][var], 0) \ for var in self.input[scen]} for scen in self.input} + # Carbon price coupling state (MSET-coupled mode only). + # Written by ftt_power.py before each solve_year(); consumed in solve_year(). + # Shape (n_reg, n_tech, 1) ÔÇö one ctax per (region, technology) in EUR2015/tCO2. + self._mset_reppx = None + + # Fuel price coupling state (MSET-coupled mode only). + # _mset_fpi is written by ftt_power.py before each solve_year() call. + # _fpix_carry accumulates the cumulative price index across years. + self._mset_fpi = None + self._fpix_carry = np.ones((len(self.titles['RTI']), len(self.titles['T2TI']), 1)) + def run(self): """ Solve model run and save results """ @@ -180,29 +217,55 @@ def solve_year(self, year, y, scenario, max_iter=1): # Run update variables, time_lags = self.update(year, y, scenario) - iter_lags = copy.deepcopy(time_lags) # Define whole period tl = self.timeline # define modules list in for possible setting.ini selection modules_list = ["FTT-P"] - # Iteration loop here - for itereration in range(max_iter): - - if "FTT-P" in self.ftt_modules: - variables = ftt_p.solve(variables, time_lags, iter_lags, - self.titles, self.conv, self.histend, tl[y], - self.domain) - - if not any(True for x in modules_list if x in self.ftt_modules): - print("Incorrect selection of modules. Check settings.ini") - - # Third, solve energy supply - # Overwrite iter_lags to be used in the next iteration round - iter_lags = copy.deepcopy(variables) -# # Print any diagnstics -# + + if "FTT-P" in self.ftt_modules: + + # Convert MINDSET ctax (EUR2015/tCO2) to CO2taxP (USD2013/tCO2). + # Use time_lags for exchange-rate snapshots: variables has them as zero (no CSV in coupled Inputs). + if self._mset_reppx is not None: + _prsc18 = time_lags.get(self._prsc_ex_var) + _ex18 = time_lags.get(self._ex_var) + if _prsc18 is not None and _ex18 is not None: + denom = (variables['PRSCX'] * _ex18 + / np.maximum(_prsc18 * variables['EXX'], 1e-10)) + variables['CO2taxP'] = (self._mset_reppx * self.carbon_price_conv_factor + / np.where(denom != 0, denom, 1.0)) + self._mset_reppx = None # consume ÔÇö must be re-set each year by ftt_power.py + + # Overwrite BCET's '22 Gamma' column with MGAM before calling FTT_Standalone. + if 'MGAM' in variables: + c2ti_m = {cat: idx for idx, cat in enumerate(self.titles['C2TI'])} + n_c2ti = len(self.titles['C2TI']) # 23 after __init__ extension + for d in (variables, time_lags): + if d['BCET'].shape[2] < n_c2ti: + ext = np.zeros((*d['BCET'].shape[:2], n_c2ti)) + ext[:, :, :d['BCET'].shape[2]] = d['BCET'] + d['BCET'] = ext + variables['BCET'][:, :, c2ti_m['22 Gamma']] = variables['MGAM'][:, :, 0] + + # FPIX is cumulative fuel price index; _fpix_carry persists it. None = standalone (no change). + fpi = self._mset_fpi if self._mset_fpi is not None else 0.0 + self._mset_fpi = None # consume ÔÇö must be re-set each year by ftt_power.py + if 'FPIX' in variables: + fpix_lag = np.where(self._fpix_carry == 0, 1.0, self._fpix_carry) + variables['FPIX'] = fpix_lag * (1.0 + fpi) + self._fpix_carry = variables['FPIX'].copy() + + # Call FTT_Standalone + variables = ftt_p_solve( + variables, time_lags, self.titles, self.histend, + tl[y], self.domain, self.power_settings + ) + + if not any(True for x in modules_list if x in self.ftt_modules): + print("Incorrect selection of modules. Check settings.ini") + return variables, time_lags def update(self, year, y, scenario): diff --git a/MINDSET_FTT_Power/SourceCode/support/input_functions.py b/MINDSET_FTT_Power/SourceCode/support/input_functions.py index 3c79d25..046793c 100644 --- a/MINDSET_FTT_Power/SourceCode/support/input_functions.py +++ b/MINDSET_FTT_Power/SourceCode/support/input_functions.py @@ -59,6 +59,12 @@ def load_data(titles, dimensions, timeline, scenarios, ftt_modules, forstart): modules_enabled = [x.strip() for x in ftt_modules.split(',')] modules_enabled += ['General'] + # Filter dims to variables whose dimensions are all resolvable in titles. + # VariableListing.csv may contain incomplete-feature variables (e.g. battery_ages, + # SCA, MAN) that reference dimension keys not yet in classification_titles.csv. + dims = {var: dimensions[var] for var in dimensions + if all(d in known_dims for d in dimensions[var])} + # Create container with the correct dimensions data = { scen : { diff --git a/MINDSET_FTT_Power/SourceCode/support/titles_functions.py b/MINDSET_FTT_Power/SourceCode/support/titles_functions.py index f1fd563..e05a318 100644 --- a/MINDSET_FTT_Power/SourceCode/support/titles_functions.py +++ b/MINDSET_FTT_Power/SourceCode/support/titles_functions.py @@ -15,43 +15,35 @@ # Third party imports -from openpyxl import load_workbook from pathlib import Path import pandas as pd def load_titles(): - # Ensure we're using consistent relative paths dir_file = os.path.dirname(os.path.realpath(__file__)) - dir_root = Path(dir_file).parents[1] - - - """ Load model classifications and titles. """ + dir_root = Path(dir_file).parents[1] - # Declare file name - titles_file = 'classification_titles.xlsx' - - # Check that classification titles workbook exists - titles_path = os.path.join(dir_root, 'Utilities', 'titles', titles_file) - if not os.path.isfile(titles_path): - print('Classification titles file not found.') + titles_path = dir_root / 'Utilities' / 'titles' / 'classification_titles.csv' + if not titles_path.is_file(): + raise FileNotFoundError(f"Classification titles file not found at: {titles_path}") - titles_wb = load_workbook(titles_path) - sheet_names = titles_wb.sheetnames - sheet_names.remove('Cover') + df = pd.read_csv(titles_path, header=None, keep_default_na=False, dtype=str) - # Iterate through worksheets and add to titles dictionary titles_dict = {} - for sheet in sheet_names: - active = titles_wb[sheet] - for column_values in active.iter_cols(min_row=1, values_only=True): - # Assigning the full names (e.g. "1 Petrol Econ") - if column_values[0] == 'Full name': # First row - titles_dict[f'{sheet}'] = column_values[1:] - # Assigning the short names (e.g. "1") - if column_values[0] == 'Short name': # First row - titles_dict[f'{sheet}_short'] = column_values[1:] + for _, row in df.iterrows(): + classification = row[0] + name_type = row[4] + values = [v for v in row.iloc[5:] if v != '' and pd.notna(v)] + cleaned = [int(v) if v.isdigit() else v for v in values] + if name_type == 'Full name': + titles_dict[classification] = tuple(cleaned) + elif name_type == 'Short name': + titles_dict[f"{classification}_short"] = tuple(cleaned) + + # '' (empty string) is used as a placeholder 4th dimension in VariableListing.csv for + # scalar/unused dims (e.g. BCET has Dim4=''). input_functions.py checks + # `all(d in known_dims for d in dims[var])` so '' must be a key in titles. + titles_dict[''] = ('',) - # Return titles dictionary return titles_dict @@ -74,6 +66,13 @@ def load_converters(): conv_dict = pd.read_excel(conv_path, sheet_name = None, index_col = 0) conv_dict.pop("Cover") + # Override T2TI_ERTI with the numbered-name mapping from ftt_t2ti_erti.csv. + # converters.xlsx uses 12 aggregate T2TI names ('Nuclear', 'Oil', ÔǪ) which diverge + # from classification_titles.csv's numbered names ('1 Nuclear', '2 Oil', ÔǪ). + # ftt_t2ti_erti.csv uses the same numbered T2TI/ERTI names as classification_titles.csv. + erti_path = os.path.join(dir_root, 'Utilities', 'ftt_t2ti_erti.csv') + if os.path.isfile(erti_path): + conv_dict['T2TI_ERTI'] = pd.read_csv(erti_path, index_col='T2TI') # Return titles dictionary return conv_dict \ No newline at end of file diff --git a/SourceCode/ftt_power.py b/SourceCode/ftt_power.py index d0e1769..a516c28 100644 --- a/SourceCode/ftt_power.py +++ b/SourceCode/ftt_power.py @@ -15,6 +15,9 @@ def solve_year_ftt(self, year, ftt_model, DYNAMIC, Scenario, ener_base, IO_model, model_start, model_end): ## FTT: Power + n_reg = len(ftt_model.titles['RTI']) + n_tech = len(ftt_model.titles['T2TI']) + _rti_short = list(ftt_model.titles['RTI_short']) # Get electricity demand and convert TJ to PJ elec_dem = ener_base.loc[ener_base.TRAD_COMM == 93].groupby('REG_imp').sum()['fuel_use'] / 1000 # Calculate FTT year index @@ -22,24 +25,26 @@ def solve_year_ftt(self, year, ftt_model, DYNAMIC, Scenario, ener_base, IO_model # Get the index of electricity elec_idx = ftt_model.titles['JTI'].index('8 Electricity') # Assign electricity demand to FTT - ftt_model.input['S0']['MEWD'][:, elec_idx, 0, y] = elec_dem[list(ftt_model.titles['RTI'])].values - # Overwrite fuel price index after 2019, when price changes are estimated in MINDSET - if year > 2019: - # Assign fuel price change to FTT - # Assess price changes by FTT fuels - # DYNAMIC['delta_price_yoy'] are domestic price changes - # The import price changes needs to be calculated from trade flows and domestic price changes - # z_bp are the monetary flows between countries - # map delta_price_yoy REG_imp to IO_model.IND_BASE REG_exp, use z_bp as weights to calculate price changes, groupby REG_imp + ftt_model.input['S0']['MEWD'][:, elec_idx, 0, y] = elec_dem[list(ftt_model.titles['RTI_short'])].values + ftt_model.input['S0']['MEWDX'][:, elec_idx, 0, y] = elec_dem[list(ftt_model.titles['RTI_short'])].values + # Overwrite fuel price index after 2019, when price changes are estimated in MINDSET. + # FPI is stored on the model object (_mset_fpi) rather than in the FTT input dict, + # so no change to FTT_Standalone's VariableListing.csv is needed. + # model_class.py converts _mset_fpi to FPIX and injects it into variables before + # calling ftt_p_solve each year. + if year > 2019: + # Assess price changes by FTT fuels. + # DYNAMIC['delta_price_yoy'] are domestic price changes. + # z_bp monetary flows weight the regional import price. fuel_pd = IO_model.IND_BASE.loc[IO_model.IND_BASE.index.get_level_values('TRAD_COMM').isin(ftt_model.ftt_fuel_converter.TRAD_COMM), 'z_bp'].copy() _dpy = self.V.read_var('delta_price_yoy', year, as_df=True).rename(columns={'delta_price_yoy': 'dp'}) exp_fuel_pd = _dpy.loc[_dpy.PROD_COMM.isin(ftt_model.ftt_tech_converter.PROD_COMM)].copy() fuel_merged = pd.merge(fuel_pd.reset_index(), exp_fuel_pd, left_on='REG_exp', right_on='REG_imp', how='inner', suffixes=('', '_y')) - + + fpi = np.zeros((n_reg, n_tech, 1)) + for tech, sectors in ftt_model.ftt_tech_converter.groupby('T2TI'): - # Get relevant sectoral data (prices and output for weighting) sec_price_chng = fuel_merged.loc[fuel_merged.PROD_COMM.isin(sectors.PROD_COMM)] - weighted_dp = (sec_price_chng.groupby("REG_imp", group_keys=False, dropna=False) .apply( lambda g: (g["dp"] * g["z_bp"]).sum() / g["z_bp"].sum() @@ -47,31 +52,33 @@ def solve_year_ftt(self, year, ftt_model, DYNAMIC, Scenario, ener_base, IO_model else 0, include_groups=False)) tech_idx = ftt_model.titles['T2TI'].index(tech) - ftt_model.input['S0']['FPI'][:, tech_idx, 0, y] = weighted_dp[list(ftt_model.titles['RTI'])].values - # Overwrite carbon prices (USD / tCO2) - for reg in ftt_model.titles['RTI']: - c_price = Scenario.carbon_tax_rate_loop(reg) - c_price = c_price.loc[c_price.PROD_COMM == 93] - for fuel, sectors in ftt_model.ftt_fuel_converter.groupby('ERTI'): - # Get relevant sectoral data (prices and output for weighting) - sec_c_price = c_price.loc[c_price.TRAD_COMM.isin(sectors.TRAD_COMM)] - if (len(sec_c_price) > 0) & (fuel in ftt_model.conv['T2TI_ERTI'].ERTI.values): - # Since carbon price is the same for all REG_exp, simply take avg. - avg_sec_c_price = sec_c_price.groupby('REG_imp')['ctax'].mean() - # Keep if weighted avreage will be needed - # weighted_dp = (sec_c_price.groupby("REG_imp", group_keys=False, dropna=False) - # .apply( - # lambda g: (g["dp"] * g["z_bp"]).sum() / g["z_bp"].sum() - # if g["z_bp"].sum() != 0 - # else 0, - # include_groups=False)) - - # Map fuel to tech - tech = ftt_model.conv['T2TI_ERTI'].index[ftt_model.conv['T2TI_ERTI'].ERTI == fuel].values[0] - tech_idx = ftt_model.titles['T2TI'].index(tech) - reg_idx = ftt_model.titles['RTI'].index(reg) - ftt_model.input['S0']['REPPX'][reg_idx, tech_idx, 0, y] = avg_sec_c_price.values[0] - + fpi[:, tech_idx, 0] = weighted_dp[_rti_short].values + + ftt_model._mset_fpi = fpi + # Carbon price: per-(region, technology) ctax in EUR2015/tCO2. + # model_class.solve_year() applies the EUR2015→USD2013 conversion using current-year FTT exchange-rate variables. + c_price_power = Scenario.tax_rate.loc[Scenario.tax_rate.PROD_COMM == 93].copy() + reppx_arr = np.zeros((n_reg, n_tech, 1)) + # converters.csv uses numbered T2TI names matching titles['T2TI'] (converters.xlsx/ftt_t2ti_erti.csv use short names — wrong for this lookup) + _conv_csv = pd.read_csv(Path('MINDSET_FTT_Power/Utilities/titles/converters.csv')) + _t2ti_list = list(ftt_model.titles['T2TI']) + _erti_to_t2ti_idxs = {} + for _, row in _conv_csv.iterrows(): + _erti_to_t2ti_idxs.setdefault(row['ERTI'], []).append(_t2ti_list.index(row['T2TI'])) + for fuel, sectors in ftt_model.ftt_fuel_converter.groupby('ERTI'): + if fuel not in _erti_to_t2ti_idxs: + continue + sec_c_price = c_price_power.loc[c_price_power.TRAD_COMM.isin(sectors.TRAD_COMM)] + if len(sec_c_price) == 0: + continue + avg_by_reg = sec_c_price.groupby('REG_imp')['ctax'].mean() + for tech_idx in _erti_to_t2ti_idxs[fuel]: + for reg_idx, reg_short in enumerate(_rti_short): + if reg_short in avg_by_reg.index: + reppx_arr[reg_idx, tech_idx, 0] = avg_by_reg[reg_short] + if np.any(reppx_arr != 0): + ftt_model._mset_reppx = reppx_arr + # else: leave _mset_reppx as None → CO2taxP stays zero (correct for no-tax years) # Solve year ftt_model.variables, ftt_model.lags = ftt_model.solve_year(year, y, ftt_model.scenarios) # Populate output container @@ -91,6 +98,9 @@ def solve_year_ftt(self, year, ftt_model, DYNAMIC, Scenario, ener_base, IO_model # Get primary energy demand values for t and t-1 ftt_energy_dem_t = ftt_model.output[ftt_model.scenarios]['MEPD'][:, :, 0, y] ftt_energy_dem_t0 = ftt_model.output[ftt_model.scenarios]['MEPD'][:, :, 0, y - 1] + # Read energy_flows once before the loop; all sector updates are applied in-place, + # then written back once after the loop (avoids N reads+writes for N fuel sectors). + _ef = self.V.read_var_df('energy_flows', year).set_index(['REG_imp', 'REG_exp', 'PROD_COMM', 'TRAD_COMM']) # Loop over supplying sectors for sector, fuels in ftt_model.ftt_fuel_converter.groupby('TRAD_COMM'): # Get indices of corresponding fuels @@ -102,7 +112,7 @@ def solve_year_ftt(self, year, ftt_model, DYNAMIC, Scenario, ener_base, IO_model ftt_energy_dem_growth = np.divide(fuel_demand_t, fuel_demand_t0, out=np.ones_like(fuel_demand_t), where=fuel_demand_t0!=0) - ftt_energy_dem_growth = pd.Series(ftt_energy_dem_growth, index = ftt_model.titles['RTI']) + ftt_energy_dem_growth = pd.Series(ftt_energy_dem_growth, index = ftt_model.titles['RTI_short']) # Filter energy data on sector power_ener_sec = power_ener_base.loc[power_ener_base.TRAD_COMM == sector].copy() power_ener_sec = power_ener_sec.rename(columns = {'fuel_use': 'fuel_use_new_adj'}) @@ -112,21 +122,23 @@ def solve_year_ftt(self, year, ftt_model, DYNAMIC, Scenario, ener_base, IO_model power_ener_sec = power_ener_sec.drop('growth', axis = 1) power_ener_sec['PROD_COMM'] = power_ener_sec['PROD_COMM'].astype(int) power_ener_sec['TRAD_COMM'] = power_ener_sec['TRAD_COMM'].astype(int) - # Make the indices aligned, apply update, write back - _ef = self.V.read_var_df('energy_flows', year).set_index(['REG_imp', 'REG_exp', 'PROD_COMM', 'TRAD_COMM']) power_ener_sec = power_ener_sec.set_index(['REG_imp', 'REG_exp', 'PROD_COMM', 'TRAD_COMM']) - # Update values based on FTT adjusted energy demand - _ef.update(power_ener_sec) # updates all overlapping columns in place - self.V.write_var_df('energy_flows', year, _ef.reset_index()) + # Accumulate updates into _ef (no write until all sectors are done) + _ef.update(power_ener_sec) + # Single write after all sector updates are applied + self.V.write_var_df('energy_flows', year, _ef.reset_index()) # Assess investment - # Calculate changes in investment - ftt_model.output[ftt_model.scenarios]['MWIY'][:, :, 0, y][:, np.newaxis, :] - ftt_model.investment[year] = (np.array(ftt_model.ftt_inv_converter[list(ftt_model.titles['T2TI'])])[np.newaxis, :, :] * - ftt_model.output[ftt_model.scenarios]['MWIY'][:, :, 0, y][:, np.newaxis, :]).sum(axis = 2) + # Reorder MWIY columns to match ftt_inv_converter, then map to MRIO sectors. + _mwiy_raw = ftt_model.output[ftt_model.scenarios]['MWIY'][:, :, 0, y] + _mwiy12 = pd.DataFrame(_mwiy_raw, columns=list(ftt_model.titles['T2TI'])) \ + .reindex(columns=list(ftt_model.ftt_inv_converter.columns)).values + ftt_model.investment[year] = ( + np.array(ftt_model.ftt_inv_converter)[np.newaxis, :, :] * + _mwiy12[:, np.newaxis, :] + ).sum(axis=2) # Convert mEUR 2010 to mUSD 2010 and then to mUSD 2019 ftt_model.investment[year] = ftt_model.investment[year] * 1.33 * 1.17 - # Export results in the last year if year == model_end: # Update scenario log diff --git a/SourceCode/government.py b/SourceCode/government.py index 40b6a19..e538af0 100644 --- a/SourceCode/government.py +++ b/SourceCode/government.py @@ -6,6 +6,7 @@ """ import pandas as pd +import numpy as np from SourceCode.utils import MRIO_df_to_vec, MRIO_vec_to_df import os @@ -181,9 +182,12 @@ def calc_gov_demand(self, government_spending_, price_index, emission_cost, new_ gov_base = gov_base[['REG_exp','REG_imp','TRAD_COMM','VIGA']].copy() consumption_shares = gov_base.merge(prices, how='left') - consumption_shares = consumption_shares.merge(emission_cost_, how='left', on=['REG_exp','REG_imp','TRAD_COMM']).fillna(0) + consumption_shares = consumption_shares.merge(emission_cost_, how='left', on=['REG_exp','REG_imp','TRAD_COMM']) + consumption_shares['emission_cost'] = consumption_shares['emission_cost'].fillna(0) + consumption_shares['price_index'] = consumption_shares['price_index'].fillna(1.0) consumption_shares['VIGA_nominal'] = consumption_shares['VIGA'] * consumption_shares['price_index'] - consumption_shares['emission_cost'] = consumption_shares['emission_cost'] / consumption_shares['VIGA_nominal'] + _denom = consumption_shares['VIGA_nominal'].replace(0, np.nan) + consumption_shares['emission_cost'] = (consumption_shares['emission_cost'] / _denom).fillna(0) consumption_shares['emission_cost'] = consumption_shares['emission_cost'].apply(lambda x: 0 if x < 0 else (3.0 if x > 3.0 else x)) consumption_shares['price_index'] = consumption_shares['price_index'] + consumption_shares['emission_cost'] consumption_shares['share'] = consumption_shares['VIGA'] / consumption_shares.groupby(['REG_imp'])['VIGA'].transform('sum') diff --git a/SourceCode/initiate_modules.py b/SourceCode/initiate_modules.py index 2ec062d..3864590 100644 --- a/SourceCode/initiate_modules.py +++ b/SourceCode/initiate_modules.py @@ -400,7 +400,7 @@ def initiate_modules(self, DYNAMIC, EXOG_VARS, MRIO_df_to_vec_DEF, MRIO_vec_to_d flow_price_impact_hh.loc[flow_price_impact_hh['flow_price_impact_hh']> 5, 'flow_price_impact_hh'] = 5.0 self.V.write_var('flow_price_impact_intermediates', year, flow_price_impact_intermediates, from_df=True) - self.V.write_var('flow_price_impact_hh', year, flow_price_impact_hh, from_df=True) + self.V.write_var('flow_price_impact_hh', year, flow_price_impact_hh.rename(columns={'PROD_COMM': 'FD'}), from_df=True) BTA_cou = BTA(Scenario, self.bta, EXOG_VARS.R, self.temp, EXOG_VARS) if year > self.model_start: @@ -409,7 +409,7 @@ def initiate_modules(self, DYNAMIC, EXOG_VARS, MRIO_df_to_vec_DEF, MRIO_vec_to_d cbam_incidence = BTA_cou.calc_cbam_incidence(EXOG_VARS.IND_BASE, DYNAMIC['carbon_content'], EXOG_VARS.HH_BASE, EXOG_VARS.FCF_BASE, EXOG_VARS.GOV_BASE) # cbam incidence ['REG_imp','REG_exp','PROD_COMM','TRAD_COMM','cbam_cost'] cbam_cost_dfs = { - 'intermediates': cbam_incidence[~cbam_incidence['PROD_COMM'].str.contains("FD")].copy(), + 'intermediates': cbam_incidence[~cbam_incidence['PROD_COMM'].str.contains("FD")].copy().astype({'PROD_COMM': int, 'TRAD_COMM': int}), 'hh': cbam_incidence[cbam_incidence['PROD_COMM']=="FD_1"].copy(), 'fcf': cbam_incidence[cbam_incidence['PROD_COMM']=="FD_4"].copy(), 'gov': cbam_incidence[cbam_incidence['PROD_COMM']=="FD_3"].copy(), @@ -672,7 +672,7 @@ def initiate_modules(self, DYNAMIC, EXOG_VARS, MRIO_df_to_vec_DEF, MRIO_vec_to_d dev_profit_rate_L1[abs(dev_profit_rate_L1)==np.inf] = 0.0 dev_profit_rate_L1[abs(dev_profit_rate_L1)>10] = 0.0 if len(dev_profit_rate_L1[abs(dev_profit_rate_L1)>10]) > 0: - print(f"profit rate issues (>10) sectors: {", ".join(map(str, np.where(dyn_qbase < 0)[0]))}") + print(f"profit rate issues (>10) sectors: {', '.join(map(str, np.where(dyn_qbase < 0)[0]))}") v_dev_profit_rate = Price_model.dp_profit_rate(dev_profit_rate_L1) dp_dev_profit_rate = Price_model.second_order_dprice(v_dev_profit_rate, year=year)['dp_full'] @@ -1033,7 +1033,8 @@ def initiate_modules(self, DYNAMIC, EXOG_VARS, MRIO_df_to_vec_DEF, MRIO_vec_to_d # Remove electricity investment from lagged investment so there is no double counting if self.ftt_run: - elec_idx = 92 * np.array(range(0, len(INV_model.R_list))) + n_sectors = len(INV_model.P_list) # 120 + elec_idx = np.arange(92, len(INV_model.R_list) * n_sectors, n_sectors) DYNAMIC['dy_inv_induced_L1'][elec_idx] = 0 dq_inv_induced, dq_inv_recyc, dq_inv_exog = IO_model.calc_dq_inv(DYNAMIC['dy_inv_induced_L1'], dy_inv_recyc, dy_inv_exog) self.V.write_var("dq_inv_exog", year, dq_inv_exog) @@ -1220,6 +1221,7 @@ def initiate_modules(self, DYNAMIC, EXOG_VARS, MRIO_df_to_vec_DEF, MRIO_vec_to_d dtax_rev, dlabor_nat = None, None A_trade_old = None cost_curves_impact_old = None + dempl_labour_supply_constraint = np.zeros_like(dempl_total) # these are the thereshold values for convergence # labor_diff, tax_diff = 0.05, 0.05 # A_trade_diff = 0.05 @@ -1231,9 +1233,9 @@ def initiate_modules(self, DYNAMIC, EXOG_VARS, MRIO_df_to_vec_DEF, MRIO_vec_to_d while self.SWITCH_WITHIN_YEAR_LOOP: iter_time = time.time() - #region 2.2.15.1_Adjust_labour_income [rgba(26,188,156,0.15)] - #region desc [rgba(26,188,156,0.50)] - # ^w [2.2.15.1] Labour income adjustment + #region 2.2.15.1_Adjust_labour_income [rgba(26,188,156,0.15)] + #region desc [rgba(26,188,156,0.50)] + # ^w [2.2.15.1] Labour income adjustment # ^w calculate delta tax revenue and take dempl_total, calculate new labour compensation based on those #endregion A_old = A_trade @@ -1250,22 +1252,25 @@ def initiate_modules(self, DYNAMIC, EXOG_VARS, MRIO_df_to_vec_DEF, MRIO_vec_to_d cost_curves_impact_old = cost_curves_impact[['REG_imp','PROD_COMM','input_cost_change']].copy() cost_curves_impact_old = cost_curves_impact_old.rename(columns={'input_cost_change':'input_cost_change_old'}) - dempl_labour_supply_constraint = np.zeros_like(dempl_total) + if iter_run > 1: + dempl_labour_supply_constraint = np.zeros_like(dempl_total) # ? calculate new emission cost based on new trade and new output Energy_emissions.update_ind_base(A_trade, self.V.read_var("output", year-1) + self.V.read_var("dq_total", year)) tax_incidence = Energy_emissions.calculate_tax_incidence() self.V.write_var_df('emission_cost_intermediates', year, tax_incidence['tax_incidence_intermediates']) - self.V.write_var_df('emission_cost_hh', year, tax_incidence['tax_incidence_hh']) - self.V.write_var_df('emission_cost_fcf', year, tax_incidence['tax_incidence_fcf']) - self.V.write_var_df('emission_cost_gov', year, tax_incidence['tax_incidence_gov']) + self.V.write_var_df('emission_cost_hh', year, tax_incidence['tax_incidence_hh'].rename(columns={'PROD_COMM': 'FD'})) + self.V.write_var_df('emission_cost_fcf', year, tax_incidence['tax_incidence_fcf'].rename(columns={'PROD_COMM': 'FD'})) + self.V.write_var_df('emission_cost_gov', year, tax_incidence['tax_incidence_gov'].rename(columns={'PROD_COMM': 'FD'})) cbam_incidence = BTA_cou.calc_cbam_incidence(ind_ener_glo, DYNAMIC['carbon_content'], EXOG_VARS.HH_BASE, EXOG_VARS.FCF_BASE, EXOG_VARS.GOV_BASE) # cbam incidence ['REG_imp','REG_exp','PROD_COMM','TRAD_COMM','cbam_cost'] - self.V.write_var_df('cbam_cost_intermediates', year, cbam_incidence[~cbam_incidence['PROD_COMM'].str.contains("FD")].copy()) - self.V.write_var_df('cbam_cost_hh', year, cbam_incidence[cbam_incidence['PROD_COMM']=="FD_1"].copy()) - self.V.write_var_df('cbam_cost_fcf', year, cbam_incidence[cbam_incidence['PROD_COMM']=="FD_4"].copy()) - self.V.write_var_df('cbam_cost_gov', year, cbam_incidence[cbam_incidence['PROD_COMM']=="FD_3"].copy()) + _cbam_interm = cbam_incidence[~cbam_incidence['PROD_COMM'].str.contains("FD")].copy() + _cbam_interm = _cbam_interm.astype({'PROD_COMM': int, 'TRAD_COMM': int}) + self.V.write_var_df('cbam_cost_intermediates', year, _cbam_interm) + self.V.write_var_df('cbam_cost_hh', year, cbam_incidence[cbam_incidence['PROD_COMM']=="FD_1"].rename(columns={'PROD_COMM':'FD'})) + self.V.write_var_df('cbam_cost_fcf', year, cbam_incidence[cbam_incidence['PROD_COMM']=="FD_4"].rename(columns={'PROD_COMM':'FD'})) + self.V.write_var_df('cbam_cost_gov', year, cbam_incidence[cbam_incidence['PROD_COMM']=="FD_3"].rename(columns={'PROD_COMM':'FD'})) # calculate within iteration delta of collected tax revenues dtax_rev = Tax_rev.calc_tax_iter_cond(emission_cost_old, self.V.read_var_df('emission_cost_intermediates', year)) @@ -1546,16 +1551,16 @@ def initiate_modules(self, DYNAMIC, EXOG_VARS, MRIO_df_to_vec_DEF, MRIO_vec_to_d dy_trade_hh = fd_trade_response['dy']['dy_trade_hh'] dy_trade_fcf = fd_trade_response['dy']['dy_trade_fcf'] dy_trade_gov = fd_trade_response['dy']['dy_trade_gov'] - + # Build new A matrix due to trade (import) substitution - + # add scenario based changes # ind_trade = IO_model.io_change(io_changes, ind_trade) A_trade = IO_model.build_A_matrix(input_df=ind_trade, variable='IO_coef_trade') IO_model.update_Leontieff(A_trade) dq_trade_eff = IO_model.calc_dq_trade((dq_tech_eff)) self.V.write_var("dq_trade_eff", year, dq_trade_eff) - + Price_model.update_A_BASE(A_trade) Price_model.calc_positive_and_negative_L() @@ -1630,7 +1635,7 @@ def initiate_modules(self, DYNAMIC, EXOG_VARS, MRIO_df_to_vec_DEF, MRIO_vec_to_d ind_ener_iter = IO_model.io_change(io_changes_, ind_trade) # ind_trade = ind_ener_iter - + # ! SET IO A_iochange = IO_model.build_A_matrix(input_df=ind_ener_iter, variable='IO_coef_trade') fd_vec_tmp = IO_model.calc_fd_vec(dq_total+IO_model.q_base_curr) @@ -1813,8 +1818,8 @@ def initiate_modules(self, DYNAMIC, EXOG_VARS, MRIO_df_to_vec_DEF, MRIO_vec_to_d print(f"--- 2.2.15 Price change relative diff.: {round(price_cond * 100, 4)}% ---") print(f"--- 2.2.15 Labour constraint diff.: {round(labor_constraint_cond * 100, 4)}% ---") print(f"--- 2.2.15 MRIO matrix relative diff.: {round(np.max(np.abs(np.nan_to_num(A_old - A_trade, nan=0))) * 100, 4)}pp ---") - - if ((labor_cond < self.COND_LABOR and tax_rev_cond < self.COND_TAX and price_cond < self.COND_PRICE and + + if ((labor_cond < self.COND_LABOR and tax_rev_cond < self.COND_TAX and price_cond < self.COND_PRICE and np.allclose(A_old, A_trade, atol=self.COND_TRADE)) or (iter_run > self.ITER_MAX)): break diff --git a/SourceCode/model_class.py b/SourceCode/model_class.py index 6ce0031..69b8773 100644 --- a/SourceCode/model_class.py +++ b/SourceCode/model_class.py @@ -81,13 +81,11 @@ from SourceCode.results import results from SourceCode.variables import variables_table from SourceCode.normal_output import normal_output_mod -from SourceCode.ftt_power import solve_year_ftt from SourceCode.cost_curves import cost_curves from SourceCode.utils import temporary_storage from SourceCode.utils import logging from SourceCode.utils import MRIO_df_to_vec, MRIO_vec_to_df, MRIO_mat_to_df from SourceCode.initiate_modules import initiate_modules -from MINDSET_FTT_Power.SourceCode.model_class import ModelRun as ftt import warnings @@ -349,7 +347,10 @@ def __init__(self): # Setup FTT:Power if self.ftt_run: - + from SourceCode.ftt_power import solve_year_ftt + from MINDSET_FTT_Power.SourceCode.model_class import ModelRun as ftt + self._solve_year_ftt = solve_year_ftt + # Instantiate the run self.ftt_model = ftt() @@ -366,23 +367,27 @@ def __init__(self): for y, year in enumerate(setup_years): # Solve year self.ftt_model.variables, self.ftt_model.lags = self.ftt_model.solve_year(year, y, self.ftt_model.scenarios) - + # Populate output container for var in self.ftt_model.variables: if 'TIME' in self.ftt_model.dims[var]: self.ftt_model.output[self.ftt_model.scenarios][var][:, :, :, y] = self.ftt_model.variables[var] else: self.ftt_model.output[self.ftt_model.scenarios][var][:, :, :, 0] = self.ftt_model.variables[var] - + # Assess investment - # Calculate changes in investment - self.ftt_model.investment[year] = (np.array(self.ftt_model.ftt_inv_converter[list(self.ftt_model.titles['T2TI'])])[np.newaxis, :, :] * - self.ftt_model.output[self.ftt_model.scenarios]['MWIY'][:, :, 0, y][:, np.newaxis, :]).sum(axis = 2) + # Reorder MWIY columns to match ftt_inv_converter, then map to MRIO sectors. + _mwiy_raw = self.ftt_model.output[self.ftt_model.scenarios]['MWIY'][:, :, 0, y] + _mwiy12 = pd.DataFrame(_mwiy_raw, columns=list(self.ftt_model.titles['T2TI'])) \ + .reindex(columns=list(self.ftt_model.ftt_inv_converter.columns)).values + self.ftt_model.investment[year] = ( + np.array(self.ftt_model.ftt_inv_converter)[np.newaxis, :, :] * + _mwiy12[:, np.newaxis, :] + ).sum(axis=2) # Convert mEUR 2010 to mUSD 2010 and then to mUSD 2019 self.ftt_model.investment[year] = self.ftt_model.investment[year] * 1.33 * 1.17 - - + #%% #region 2_Solving_the_model [rgba(52,152,219,0.10)] @@ -425,7 +430,7 @@ def __init__(self): 'dq_exog_gov_prev': np.zeros((len(self.EXOG_VARS.R)*len(self.EXOG_VARS.P))), 'dq_supply_constraint_no_empl': np.zeros((len(self.EXOG_VARS.R)*len(self.EXOG_VARS.P))), 'fuel_price': pd.DataFrame(), - 'cbam_incidence': {} + 'cbam_incidence': {}, } self.DYNAMIC['PROJECTION_OUTPUT'] = self.EXOG_VARS.PROJECTION_OUTPUT @@ -1064,13 +1069,9 @@ def solve_year(self, year, DYNAMIC, EXOG_VARS, CALIBRATING, MRIO_df_to_vec_DEF, if self.ftt_run: ## FTT: Power - ftt_model, DYNAMIC = solve_year_ftt(self, year, ftt_model, DYNAMIC, Scenario, - ener_base, IO_model, self.model_start, self.model_end) + ftt_model, DYNAMIC = self._solve_year_ftt(self, year, ftt_model, DYNAMIC, Scenario, + ener_base, IO_model, self.model_start, self.model_end) - # for i in range(len(ftt_model.titles['RTI'])): - # print(ftt_model.titles['RTI'][i]) - # pd.DataFrame(ftt_model.output['S0']['MEWG'][i, :, 0, :], index = ftt_model.titles['T2TI'], columns = range(2010, 2051)).T.plot() - ## EMISSIONS Energy_emissions = ener_balance(EXOG_VARS, Scenario, self.refining_sectors)