# coding: utf-8 import logging logger = logging.getLogger(__name__) import pandas as pd idx = pd.IndexSlice import numpy as np import scipy as sp import xarray as xr import re, os from six import iteritems, string_types import pypsa import yaml import pytz from vresutils.costdata import annuity from scipy.stats import beta from build_energy_totals import build_eea_co2, build_eurostat_co2, build_co2_totals #First tell PyPSA that links can have multiple outputs by #overriding the component_attrs. This can be done for #as many buses as you need with format busi for i = 2,3,4,5,.... #See https://pypsa.org/doc/components.html#link-with-multiple-outputs-or-inputs override_component_attrs = pypsa.descriptors.Dict({k : v.copy() for k,v in pypsa.components.component_attrs.items()}) override_component_attrs["Link"].loc["bus2"] = ["string",np.nan,np.nan,"2nd bus","Input (optional)"] override_component_attrs["Link"].loc["bus3"] = ["string",np.nan,np.nan,"3rd bus","Input (optional)"] override_component_attrs["Link"].loc["bus4"] = ["string",np.nan,np.nan,"4th bus","Input (optional)"] override_component_attrs["Link"].loc["efficiency2"] = ["static or series","per unit",1.,"2nd bus efficiency","Input (optional)"] override_component_attrs["Link"].loc["efficiency3"] = ["static or series","per unit",1.,"3rd bus efficiency","Input (optional)"] override_component_attrs["Link"].loc["efficiency4"] = ["static or series","per unit",1.,"4th bus efficiency","Input (optional)"] override_component_attrs["Link"].loc["p2"] = ["series","MW",0.,"2nd bus output","Output"] override_component_attrs["Link"].loc["p3"] = ["series","MW",0.,"3rd bus output","Output"] override_component_attrs["Link"].loc["p4"] = ["series","MW",0.,"4th bus output","Output"] override_component_attrs["Link"].loc["build_year"] = ["integer","year",np.nan,"build year","Input (optional)"] override_component_attrs["Link"].loc["lifetime"] = ["float","years",np.nan,"lifetime","Input (optional)"] override_component_attrs["Generator"].loc["build_year"] = ["integer","year",np.nan,"build year","Input (optional)"] override_component_attrs["Generator"].loc["lifetime"] = ["float","years",np.nan,"lifetime","Input (optional)"] override_component_attrs["Store"].loc["build_year"] = ["integer","year",np.nan,"build year","Input (optional)"] override_component_attrs["Store"].loc["lifetime"] = ["float","years",np.nan,"lifetime","Input (optional)"] def co2_emissions_year(cts, opts, year): """ calculate co2 emissions in one specific year (e.g. 1990 or 2018). """ eea_co2 = build_eea_co2(year) #TODO: read Eurostat data from year>2014, this only affects the estimation of # CO2 emissions for "BA","RS","AL","ME","MK" if year > 2014: eurostat_co2 = build_eurostat_co2(year=2014) else: eurostat_co2 = build_eurostat_co2(year) co2_totals=build_co2_totals(eea_co2, eurostat_co2, year) co2_emissions = co2_totals.loc[cts, "electricity"].sum() if "T" in opts: co2_emissions += co2_totals.loc[cts, [i+ " non-elec" for i in ["rail","road"]]].sum().sum() if "H" in opts: co2_emissions += co2_totals.loc[cts, [i+ " non-elec" for i in ["residential","services"]]].sum().sum() if "I" in opts: co2_emissions += co2_totals.loc[cts, ["industrial non-elec","industrial processes", "domestic aviation","international aviation", "domestic navigation","international navigation"]].sum().sum() co2_emissions *=0.001 #MtCO2 to GtCO2 return co2_emissions def build_carbon_budget(o): #distribute carbon budget following beta or exponential transition path if "be" in o: #beta decay carbon_budget = float(o[o.find("cb")+2:o.find("be")]) be=float(o[o.find("be")+2:]) if "ex" in o: #exponential decay carbon_budget = float(o[o.find("cb")+2:o.find("ex")]) r=float(o[o.find("ex")+2:]) pop_layout = pd.read_csv(snakemake.input.clustered_pop_layout, index_col=0) pop_layout["ct"] = pop_layout.index.str[:2] cts = pop_layout.ct.value_counts().index e_1990 = co2_emissions_year(cts, opts, year=1990) #emissions at the beginning of the path (last year available 2018) e_0 = co2_emissions_year(cts, opts, year=2018) #emissions in 2019 and 2020 assumed equal to 2018 and substracted carbon_budget -= 2*e_0 planning_horizons = snakemake.config['scenario']['planning_horizons'] CO2_CAP = pd.DataFrame(index = pd.Series(data=planning_horizons, name='planning_horizon'), columns=pd.Series(data=[], name='paths', dtype='float')) t_0 = planning_horizons[0] if "be" in o: #beta decay t_f = t_0 + (2*carbon_budget/e_0).round(0) # final year in the path #emissions (relative to 1990) CO2_CAP[o] = [(e_0/e_1990)*(1-beta.cdf((t-t_0)/(t_f-t_0), be, be)) for t in planning_horizons] if "ex" in o: #exponential decay without delay T=carbon_budget/e_0 m=(1+np.sqrt(1+r*T))/T CO2_CAP[o] = [(e_0/e_1990)*(1+(m+r)*(t-t_0))*np.exp(-m*(t-t_0)) for t in planning_horizons] CO2_CAP.to_csv(path_cb + 'carbon_budget_distribution.csv', sep=',', line_terminator='\n', float_format='%.3f') countries=pd.Series(data=cts) countries.to_csv(path_cb + 'countries.csv', sep=',', line_terminator='\n', float_format='%.3f') def add_lifetime_wind_solar(n): """ Add lifetime for solar and wind generators """ for carrier in ['solar', 'onwind', 'offwind-dc', 'offwind-ac']: carrier_name='offwind' if carrier in ['offwind-dc', 'offwind-ac'] else carrier n.generators.loc[[index for index in n.generators.index.to_list() if carrier in index], 'lifetime']=costs.at[carrier_name,'lifetime'] def update_wind_solar_costs(n,costs): """ Update costs for wind and solar generators added with pypsa-eur to those cost in the planning year """ #NB: solar costs are also manipulated for rooftop #when distribution grid is inserted n.generators.loc[n.generators.carrier=='solar','capital_cost'] = costs.at['solar-utility', 'fixed'] n.generators.loc[n.generators.carrier=='onwind','capital_cost'] = costs.at['onwind', 'fixed'] #for offshore wind, need to calculated connection costs #assign clustered bus #map initial network -> simplified network busmap_s = pd.read_csv(snakemake.input.busmap_s, index_col=0).squeeze() busmap_s.index = busmap_s.index.astype(str) busmap_s = busmap_s.astype(str) #map simplified network -> clustered network busmap = pd.read_csv(snakemake.input.busmap, index_col=0).squeeze() busmap.index = busmap.index.astype(str) busmap = busmap.astype(str) #map initial network -> clustered network clustermaps = busmap_s.map(busmap) #code adapted from pypsa-eur/scripts/add_electricity.py for connection in ['dc','ac']: tech = "offwind-" + connection profile = snakemake.input['profile_offwind_' + connection] with xr.open_dataset(profile) as ds: underwater_fraction = ds['underwater_fraction'].to_pandas() connection_cost = (snakemake.config['costs']['lines']['length_factor'] * ds['average_distance'].to_pandas() * (underwater_fraction * costs.at[tech + '-connection-submarine', 'fixed'] + (1. - underwater_fraction) * costs.at[tech + '-connection-underground', 'fixed'])) #convert to aggregated clusters with weighting weight = ds['weight'].to_pandas() #e.g. clusters == 37m means that VRE generators are left #at clustering of simplified network, but that they are #connected to 37-node network if snakemake.wildcards.clusters[-1:] == "m": genmap = busmap_s else: genmap = clustermaps connection_cost = (connection_cost*weight).groupby(genmap).sum()/weight.groupby(genmap).sum() capital_cost = (costs.at['offwind', 'fixed'] + costs.at[tech + '-station', 'fixed'] + connection_cost) logger.info("Added connection cost of {:0.0f}-{:0.0f} Eur/MW/a to {}" .format(connection_cost[0].min(), connection_cost[0].max(), tech)) n.generators.loc[n.generators.carrier==tech,'capital_cost'] = capital_cost.rename(index=lambda node: node + ' ' + tech) def add_carrier_buses(n, carriers): """ Add buses to connect e.g. coal, nuclear and oil plants """ for carrier in carriers: n.add("Carrier", carrier) #use madd to get location inserted n.madd("Bus", ["EU " + carrier], location="EU", carrier=carrier) #use madd to get carrier inserted n.madd("Store", ["EU " + carrier + " Store"], bus=["EU " + carrier], e_nom_extendable=True, e_cyclic=True, carrier=carrier, capital_cost=0.) #could correct to e.g. 0.2 EUR/kWh * annuity and O&M n.add("Generator", "EU " + carrier, bus="EU " + carrier, p_nom_extendable=True, carrier=carrier, capital_cost=0., marginal_cost=costs.at[carrier,'fuel']) def remove_elec_base_techs(n): """remove conventional generators (e.g. OCGT) and storage units (e.g. batteries and H2) from base electricity-only network, since they're added here differently using links """ for c in n.iterate_components(snakemake.config["pypsa_eur"]): to_keep = snakemake.config["pypsa_eur"][c.name] to_remove = pd.Index(c.df.carrier.unique())^to_keep print("Removing",c.list_name,"with carrier",to_remove) names = c.df.index[c.df.carrier.isin(to_remove)] print(names) n.mremove(c.name, names) n.carriers.drop(to_remove, inplace=True, errors="ignore") def add_co2_tracking(n): #minus sign because opposite to how fossil fuels used: #CH4 burning puts CH4 down, atmosphere up n.add("Carrier","co2", co2_emissions=-1.) #this tracks CO2 in the atmosphere n.madd("Bus", ["co2 atmosphere"], location="EU", carrier="co2") #NB: can also be negative n.madd("Store",["co2 atmosphere"], e_nom_extendable=True, e_min_pu=-1, carrier="co2", bus="co2 atmosphere") #this tracks CO2 stored, e.g. underground n.madd("Bus", ["co2 stored"], location="EU", carrier="co2 stored") n.madd("Store",["co2 stored"], e_nom_extendable=True, e_nom_max=options['co2_sequestration_potential']*1e6, capital_cost=options['co2_sequestration_cost'], carrier="co2 stored", bus="co2 stored") if options['co2_vent']: n.madd("Link",["co2 vent"], bus0="co2 stored", bus1="co2 atmosphere", carrier="co2 vent", efficiency=1., p_nom_extendable=True) def add_dac(n): heat_buses = n.buses.index[n.buses.carrier.isin(["urban central heat", "services urban decentral heat"])] locations = n.buses.location[heat_buses] n.madd("Link", locations, suffix=" DAC", bus0="co2 atmosphere", bus1="co2 stored", bus2=locations.values, bus3=heat_buses, carrier="DAC", capital_cost=costs.at['direct air capture','fixed'], efficiency=1., efficiency2=-(costs.at['direct air capture','electricity-input'] + costs.at['direct air capture','compression-electricity-input']), efficiency3=-(costs.at['direct air capture','heat-input'] - costs.at['direct air capture','compression-heat-output']), p_nom_extendable=True, lifetime=costs.at['direct air capture','lifetime']) def add_co2limit(n, Nyears=1.,limit=0.): cts = pop_layout.ct.value_counts().index co2_limit = co2_totals.loc[cts, "electricity"].sum() if "T" in opts: co2_limit += co2_totals.loc[cts, [i+ " non-elec" for i in ["rail","road"]]].sum().sum() if "H" in opts: co2_limit += co2_totals.loc[cts, [i+ " non-elec" for i in ["residential","services"]]].sum().sum() if "I" in opts: co2_limit += co2_totals.loc[cts, ["industrial non-elec","industrial processes", "domestic aviation","international aviation", "domestic navigation","international navigation"]].sum().sum() co2_limit *= limit*Nyears n.add("GlobalConstraint", "CO2Limit", carrier_attribute="co2_emissions", sense="<=", constant=co2_limit) def add_emission_prices(n, emission_prices=None, exclude_co2=False): assert False, "Needs to be fixed, adds NAN" if emission_prices is None: emission_prices = snakemake.config['costs']['emission_prices'] if exclude_co2: emission_prices.pop('co2') ep = (pd.Series(emission_prices).rename(lambda x: x+'_emissions') * n.carriers).sum(axis=1) n.generators['marginal_cost'] += n.generators.carrier.map(ep) n.storage_units['marginal_cost'] += n.storage_units.carrier.map(ep) def set_line_s_max_pu(n): # set n-1 security margin to 0.5 for 37 clusters and to 0.7 from 200 clusters # 128 reproduces 98% of line volume in TWkm, but clustering distortions inside node n_clusters = len(n.buses.index[n.buses.carrier == "AC"]) s_max_pu = np.clip(0.5 + 0.2 * (n_clusters - 37) / (200 - 37), 0.5, 0.7) n.lines['s_max_pu'] = s_max_pu dc_b = n.links.carrier == 'DC' n.links.loc[dc_b, 'p_max_pu'] = snakemake.config['links']['p_max_pu'] n.links.loc[dc_b, 'p_min_pu'] = - snakemake.config['links']['p_max_pu'] def set_line_volume_limit(n, lv): dc_b = n.links.carrier == 'DC' if lv != "opt": lv = float(lv) # Either line_volume cap or cost n.lines['capital_cost'] = 0. n.links.loc[dc_b,'capital_cost'] = 0. else: n.lines['capital_cost'] = (n.lines['length'] * costs.at['HVAC overhead', 'fixed']) #add HVDC inverter post factor, to maintain consistency with LV limit n.links.loc[dc_b, 'capital_cost'] = (n.links.loc[dc_b, 'length'] * costs.at['HVDC overhead', 'fixed'])# + #costs.at['HVDC inverter pair', 'fixed']) if lv != 1.0: lines_s_nom = n.lines.s_nom.where( n.lines.type == '', np.sqrt(3) * n.lines.num_parallel * n.lines.type.map(n.line_types.i_nom) * n.lines.bus0.map(n.buses.v_nom) ) n.lines['s_nom_min'] = lines_s_nom n.links.loc[dc_b,'p_nom_min'] = n.links['p_nom'] n.lines['s_nom_extendable'] = True n.links.loc[dc_b,'p_nom_extendable'] = True if lv != "opt": n.line_volume_limit = lv * ((lines_s_nom * n.lines['length']).sum() + n.links.loc[dc_b].eval('p_nom * length').sum()) return n def average_every_nhours(n, offset): logger.info('Resampling the network to {}'.format(offset)) m = n.copy(with_time=False) #fix copying of network attributes #copied from pypsa/io.py, should be in pypsa/components.py#Network.copy() allowed_types = (float,int,bool) + string_types + tuple(np.typeDict.values()) attrs = dict((attr, getattr(n, attr)) for attr in dir(n) if (not attr.startswith("__") and isinstance(getattr(n,attr), allowed_types))) for k,v in iteritems(attrs): setattr(m,k,v) snapshot_weightings = n.snapshot_weightings.resample(offset).sum() m.set_snapshots(snapshot_weightings.index) m.snapshot_weightings = snapshot_weightings for c in n.iterate_components(): pnl = getattr(m, c.list_name+"_t") for k, df in iteritems(c.pnl): if not df.empty: if c.list_name == "stores" and k == "e_max_pu": pnl[k] = df.resample(offset).min() elif c.list_name == "stores" and k == "e_min_pu": pnl[k] = df.resample(offset).max() else: pnl[k] = df.resample(offset).mean() return m def generate_periodic_profiles(dt_index=pd.date_range("2011-01-01 00:00","2011-12-31 23:00",freq="H",tz="UTC"), nodes=[], weekly_profile=range(24*7)): """Give a 24*7 long list of weekly hourly profiles, generate this for each country for the period dt_index, taking account of time zones and Summer Time. """ weekly_profile = pd.Series(weekly_profile,range(24*7)) week_df = pd.DataFrame(index=dt_index,columns=nodes) for ct in nodes: week_df[ct] = [24*dt.weekday()+dt.hour for dt in dt_index.tz_convert(pytz.timezone(timezone_mappings[ct[:2]]))] week_df[ct] = week_df[ct].map(weekly_profile) return week_df def shift_df(df,hours=1): """Works both on Series and DataFrame""" df = df.copy() df.values[:] = np.concatenate([df.values[-hours:], df.values[:-hours]]) return df def transport_degree_factor(temperature,deadband_lower=15,deadband_upper=20, lower_degree_factor=0.5, upper_degree_factor=1.6): """Work out how much energy demand in vehicles increases due to heating and cooling. There is a deadband where there is no increase. Degree factors are % increase in demand compared to no heating/cooling fuel consumption. Returns per unit increase in demand for each place and time """ dd = temperature.copy() dd[(temperature > deadband_lower) & (temperature < deadband_upper)] = 0. dd[temperature < deadband_lower] = lower_degree_factor/100.*(deadband_lower-temperature[temperature < deadband_lower]) dd[temperature > deadband_upper] = upper_degree_factor/100.*(temperature[temperature > deadband_upper]-deadband_upper) return dd def prepare_data(network): ############## #Heating ############## ashp_cop = xr.open_dataarray(snakemake.input.cop_air_total).T.to_pandas().reindex(index=network.snapshots) gshp_cop = xr.open_dataarray(snakemake.input.cop_soil_total).T.to_pandas().reindex(index=network.snapshots) solar_thermal = xr.open_dataarray(snakemake.input.solar_thermal_total).T.to_pandas().reindex(index=network.snapshots) #1e3 converts from W/m^2 to MW/(1000m^2) = kW/m^2 solar_thermal = options['solar_cf_correction'] * solar_thermal/1e3 energy_totals = pd.read_csv(snakemake.input.energy_totals_name,index_col=0) nodal_energy_totals = energy_totals.loc[pop_layout.ct].fillna(0.) nodal_energy_totals.index = pop_layout.index nodal_energy_totals = nodal_energy_totals.multiply(pop_layout.fraction,axis=0) #copy forward the daily average heat demand into each hour, so it can be multipled by the intraday profile daily_space_heat_demand = xr.open_dataarray(snakemake.input.heat_demand_total).T.to_pandas().reindex(index=network.snapshots, method="ffill") intraday_profiles = pd.read_csv(snakemake.input.heat_profile,index_col=0) sectors = ["residential","services"] uses = ["water","space"] heat_demand = {} electric_heat_supply = {} for sector in sectors: for use in uses: intraday_year_profile = generate_periodic_profiles(daily_space_heat_demand.index.tz_localize("UTC"), nodes=daily_space_heat_demand.columns, weekly_profile=(list(intraday_profiles["{} {} weekday".format(sector,use)])*5 + list(intraday_profiles["{} {} weekend".format(sector,use)])*2)).tz_localize(None) if use == "space": heat_demand_shape = daily_space_heat_demand*intraday_year_profile else: heat_demand_shape = intraday_year_profile heat_demand["{} {}".format(sector,use)] = (heat_demand_shape/heat_demand_shape.sum()).multiply(nodal_energy_totals["total {} {}".format(sector,use)])*1e6 electric_heat_supply["{} {}".format(sector,use)] = (heat_demand_shape/heat_demand_shape.sum()).multiply(nodal_energy_totals["electricity {} {}".format(sector,use)])*1e6 heat_demand = pd.concat(heat_demand,axis=1) electric_heat_supply = pd.concat(electric_heat_supply,axis=1) #subtract from electricity load since heat demand already in heat_demand electric_nodes = n.loads.index[n.loads.carrier == "electricity"] n.loads_t.p_set[electric_nodes] = n.loads_t.p_set[electric_nodes] - electric_heat_supply.groupby(level=1,axis=1).sum()[electric_nodes] ############## #Transport ############## ## Get overall demand curve for all vehicles traffic = pd.read_csv(os.path.join(snakemake.input.traffic_data,"KFZ__count"), skiprows=2)["count"] #Generate profiles transport_shape = generate_periodic_profiles(dt_index=network.snapshots.tz_localize("UTC"), nodes=pop_layout.index, weekly_profile=traffic.values).tz_localize(None) transport_shape = transport_shape/transport_shape.sum() transport_data = pd.read_csv(snakemake.input.transport_name, index_col=0) nodal_transport_data = transport_data.loc[pop_layout.ct].fillna(0.) nodal_transport_data.index = pop_layout.index nodal_transport_data["number cars"] = pop_layout["fraction"]*nodal_transport_data["number cars"] nodal_transport_data.loc[nodal_transport_data["average fuel efficiency"] == 0.,"average fuel efficiency"] = transport_data["average fuel efficiency"].mean() #electric motors are more efficient, so alter transport demand #kWh/km from EPA https://www.fueleconomy.gov/feg/ for Tesla Model S plug_to_wheels_eta = 0.20 battery_to_wheels_eta = plug_to_wheels_eta*0.9 efficiency_gain = nodal_transport_data["average fuel efficiency"]/battery_to_wheels_eta #get heating demand for correction to demand time series temperature = xr.open_dataarray(snakemake.input.temp_air_total).T.to_pandas() #correction factors for vehicle heating dd_ICE = transport_degree_factor(temperature, options['transport_heating_deadband_lower'], options['transport_heating_deadband_upper'], options['ICE_lower_degree_factor'], options['ICE_upper_degree_factor']) dd_EV = transport_degree_factor(temperature, options['transport_heating_deadband_lower'], options['transport_heating_deadband_upper'], options['EV_lower_degree_factor'], options['EV_upper_degree_factor']) #divide out the heating/cooling demand from ICE totals ICE_correction = (transport_shape*(1+dd_ICE)).sum()/transport_shape.sum() transport = (transport_shape.multiply(nodal_energy_totals["total road"] + nodal_energy_totals["total rail"] - nodal_energy_totals["electricity rail"])*1e6*Nyears).divide(efficiency_gain*ICE_correction) #multiply back in the heating/cooling demand for EVs transport = transport.multiply(1+dd_EV) ## derive plugged-in availability for PKW's (cars) traffic = pd.read_csv(os.path.join(snakemake.input.traffic_data,"Pkw__count"), skiprows=2)["count"] avail_max = 0.95 avail_mean = 0.8 avail = avail_max - (avail_max - avail_mean)*(traffic - traffic.min())/(traffic.mean() - traffic.min()) avail_profile = generate_periodic_profiles(dt_index=network.snapshots.tz_localize("UTC"), nodes=pop_layout.index, weekly_profile=avail.values).tz_localize(None) dsm_week = np.zeros((24*7,)) dsm_week[(np.arange(0,7,1)*24+options['bev_dsm_restriction_time'])] = options['bev_dsm_restriction_value'] dsm_profile = generate_periodic_profiles(dt_index=network.snapshots.tz_localize("UTC"), nodes=pop_layout.index, weekly_profile=dsm_week).tz_localize(None) ############### #CO2 ############### #1e6 to convert Mt to tCO2 co2_totals = 1e6*pd.read_csv(snakemake.input.co2_totals_name,index_col=0) return nodal_energy_totals, heat_demand, ashp_cop, gshp_cop, solar_thermal, transport, avail_profile, dsm_profile, co2_totals, nodal_transport_data def prepare_costs(cost_file, USD_to_EUR, discount_rate, Nyears, lifetime): #set all asset costs and other parameters costs = pd.read_csv(cost_file,index_col=list(range(2))).sort_index() #correct units to MW and EUR costs.loc[costs.unit.str.contains("/kW"),"value"]*=1e3 costs.loc[costs.unit.str.contains("USD"),"value"]*=USD_to_EUR #min_count=1 is important to generate NaNs which are then filled by fillna costs = costs.loc[:, "value"].unstack(level=1).groupby("technology").sum(min_count=1) costs = costs.fillna({"CO2 intensity" : 0, "FOM" : 0, "VOM" : 0, "discount rate" : discount_rate, "efficiency" : 1, "fuel" : 0, "investment" : 0, "lifetime" : lifetime }) costs["fixed"] = [(annuity(v["lifetime"],v["discount rate"])+v["FOM"]/100.)*v["investment"]*Nyears for i,v in costs.iterrows()] return costs def add_generation(network): print("adding electricity generation") nodes = pop_layout.index conventionals = [("OCGT","gas")] for generator,carrier in [("OCGT","gas")]: network.add("Carrier", carrier) network.madd("Bus", ["EU " + carrier], location="EU", carrier=carrier) #use madd to get carrier inserted network.madd("Store", ["EU " + carrier + " Store"], bus=["EU " + carrier], e_nom_extendable=True, e_cyclic=True, carrier=carrier, capital_cost=0.) #could correct to e.g. 0.2 EUR/kWh * annuity and O&M network.add("Generator", "EU " + carrier, bus="EU " + carrier, p_nom_extendable=True, carrier=carrier, capital_cost=0., marginal_cost=costs.at[carrier,'fuel']) network.madd("Link", nodes + " " + generator, bus0=["EU " + carrier]*len(nodes), bus1=nodes, bus2="co2 atmosphere", marginal_cost=costs.at[generator,'efficiency']*costs.at[generator,'VOM'], #NB: VOM is per MWel capital_cost=costs.at[generator,'efficiency']*costs.at[generator,'fixed'], #NB: fixed cost is per MWel p_nom_extendable=True, carrier=generator, efficiency=costs.at[generator,'efficiency'], efficiency2=costs.at[carrier,'CO2 intensity'], lifetime=costs.at[generator,'lifetime']) def add_wave(network, wave_cost_factor): wave_fn = "data/WindWaveWEC_GLTB.xlsx" locations = ["FirthForth","Hebrides"] #in kW capacity = pd.Series([750,1000,600],["Attenuator","F2HB","MultiPA"]) #in EUR/MW costs = wave_cost_factor*pd.Series([2.5,2,1.5],["Attenuator","F2HB","MultiPA"])*1e6 sheets = {} for l in locations: sheets[l] = pd.read_excel(wave_fn, index_col=0,skiprows=[0],parse_dates=True, sheet_name=l) to_drop = ["Vestas 3MW","Vestas 8MW"] wave = pd.concat([sheets[l].drop(to_drop,axis=1).divide(capacity,axis=1) for l in locations], keys=locations, axis=1) for wave_type in costs.index: n.add("Generator", "Hebrides "+wave_type, bus="GB4 0", p_nom_extendable=True, carrier="wave", capital_cost=(annuity(25,0.07)+0.03)*costs[wave_type], p_max_pu=wave["Hebrides",wave_type]) def insert_electricity_distribution_grid(network): print("Inserting electricity distribution grid with investment cost factor of", snakemake.config["sector"]['electricity_distribution_grid_cost_factor']) nodes = pop_layout.index network.madd("Bus", nodes+ " low voltage", location=nodes, carrier="low voltage") network.madd("Link", nodes + " electricity distribution grid", bus0=nodes, bus1=nodes + " low voltage", p_nom_extendable=True, p_min_pu=-1, carrier="electricity distribution grid", efficiency=1, marginal_cost=0, lifetime=costs.at['electricity distribution grid','lifetime'], capital_cost=costs.at['electricity distribution grid','fixed']*snakemake.config["sector"]['electricity_distribution_grid_cost_factor']) #this catches regular electricity load and "industry electricity" loads = network.loads.index[network.loads.carrier.str.contains("electricity")] network.loads.loc[loads,"bus"] += " low voltage" bevs = network.links.index[network.links.carrier == "BEV charger"] network.links.loc[bevs,"bus0"] += " low voltage" v2gs = network.links.index[network.links.carrier == "V2G"] network.links.loc[v2gs,"bus1"] += " low voltage" hps = network.links.index[network.links.carrier.str.contains("heat pump")] network.links.loc[hps,"bus0"] += " low voltage" rh = network.links.index[network.links.carrier.str.contains("resistive heater")] network.links.loc[rh, "bus0"] += " low voltage" mchp = network.links.index[network.links.carrier.str.contains("micro gas")] network.links.loc[mchp, "bus1"] += " low voltage" #set existing solar to cost of utility cost rather the 50-50 rooftop-utility solar = network.generators.index[network.generators.carrier == "solar"] network.generators.loc[solar, "capital_cost"] = costs.at['solar-utility', 'fixed'] if snakemake.wildcards.clusters[-1:] == "m": pop_solar = simplified_pop_layout.total.rename(index = lambda x: x + " solar") else: pop_solar = pop_layout.total.rename(index = lambda x: x + " solar") # add max solar rooftop potential assuming 0.1 kW/m2 and 10 m2/person, #i.e. 1 kW/person (population data is in thousands of people) so we get MW potential = 0.1*10*pop_solar network.madd("Generator", solar, suffix=" rooftop", bus=network.generators.loc[solar, "bus"] + " low voltage", carrier="solar rooftop", p_nom_extendable=True, p_nom_max=potential, marginal_cost=network.generators.loc[solar, 'marginal_cost'], capital_cost=costs.at['solar-rooftop', 'fixed'], efficiency=network.generators.loc[solar, 'efficiency'], p_max_pu=network.generators_t.p_max_pu[solar]) network.add("Carrier","home battery") network.madd("Bus", nodes + " home battery", location=nodes, carrier="home battery") network.madd("Store", nodes + " home battery", bus=nodes + " home battery", e_cyclic=True, e_nom_extendable=True, carrier="home battery", capital_cost=costs.at['battery storage','fixed'], lifetime=costs.at['battery storage','lifetime']) network.madd("Link", nodes + " home battery charger", bus0=nodes + " low voltage", bus1=nodes + " home battery", carrier="home battery charger", efficiency=costs.at['battery inverter','efficiency']**0.5, capital_cost=costs.at['battery inverter','fixed'], p_nom_extendable=True, lifetime=costs.at['battery inverter','lifetime']) network.madd("Link", nodes + " home battery discharger", bus0=nodes + " home battery", bus1=nodes + " low voltage", carrier="home battery discharger", efficiency=costs.at['battery inverter','efficiency']**0.5, marginal_cost=options['marginal_cost_storage'], p_nom_extendable=True, lifetime=costs.at['battery inverter','lifetime']) def insert_gas_distribution_costs(network): f_costs = options['gas_distribution_grid_cost_factor'] print("Inserting gas distribution grid with investment cost\ factor of", f_costs) # gas boilers gas_b = network.links[network.links.carrier.str.contains("gas boiler") & (~network.links.carrier.str.contains("urban central"))].index network.links.loc[gas_b, "capital_cost"] += costs.loc['electricity distribution grid']["fixed"] * f_costs # micro CHPs mchp = network.links.index[network.links.carrier.str.contains("micro gas")] network.links.loc[mchp, "capital_cost"] += costs.loc['electricity distribution grid']["fixed"] * f_costs def add_electricity_grid_connection(network): carriers = ["onwind","solar"] gens = network.generators.index[network.generators.carrier.isin(carriers)] network.generators.loc[gens,"capital_cost"] += costs.at['electricity grid connection','fixed'] def add_storage(network): print("adding electricity storage") nodes = pop_layout.index network.add("Carrier","H2") network.madd("Bus", nodes+ " H2", location=nodes, carrier="H2") network.madd("Link", nodes + " H2 Electrolysis", bus1=nodes + " H2", bus0=nodes, p_nom_extendable=True, carrier="H2 Electrolysis", efficiency=costs.at["electrolysis","efficiency"], capital_cost=costs.at["electrolysis","fixed"], lifetime=costs.at['electrolysis','lifetime']) network.madd("Link", nodes + " H2 Fuel Cell", bus0=nodes + " H2", bus1=nodes, p_nom_extendable=True, carrier ="H2 Fuel Cell", efficiency=costs.at["fuel cell","efficiency"], capital_cost=costs.at["fuel cell","fixed"]*costs.at["fuel cell","efficiency"], #NB: fixed cost is per MWel lifetime=costs.at['fuel cell','lifetime']) cavern_nodes = pd.DataFrame() if options['hydrogen_underground_storage']: h2_salt_cavern_potential = pd.read_csv(snakemake.input.h2_cavern, index_col=0,squeeze=True) h2_cavern_ct = h2_salt_cavern_potential[~h2_salt_cavern_potential.isna()] cavern_nodes = pop_layout[pop_layout.ct.isin(h2_cavern_ct.index)] h2_capital_cost = costs.at["hydrogen storage underground", "fixed"] # assumptions: weight storage potential in a country by population # TODO: fix with real geographic potentials #convert TWh to MWh with 1e6 h2_pot = h2_cavern_ct.loc[cavern_nodes.ct] h2_pot.index = cavern_nodes.index h2_pot = h2_pot * cavern_nodes.fraction * 1e6 network.madd("Store", cavern_nodes.index + " H2 Store", bus=cavern_nodes.index + " H2", e_nom_extendable=True, e_nom_max=h2_pot.values, e_cyclic=True, carrier="H2 Store", capital_cost=h2_capital_cost) # hydrogen stored overground h2_capital_cost = costs.at["hydrogen storage tank", "fixed"] nodes_overground = nodes ^ cavern_nodes.index network.madd("Store", nodes_overground + " H2 Store", bus=nodes_overground + " H2", e_nom_extendable=True, e_cyclic=True, carrier="H2 Store", capital_cost=h2_capital_cost) h2_links = pd.DataFrame(columns=["bus0","bus1","length"]) prefix = "H2 pipeline " connector = " -> " attrs = ["bus0","bus1","length"] candidates = pd.concat([network.lines[attrs],network.links.loc[network.links.carrier == "DC",attrs]], keys=["lines","links"]) for candidate in candidates.index: buses = [candidates.at[candidate,"bus0"],candidates.at[candidate,"bus1"]] buses.sort() name = prefix + buses[0] + connector + buses[1] if name not in h2_links.index: h2_links.at[name,"bus0"] = buses[0] h2_links.at[name,"bus1"] = buses[1] h2_links.at[name,"length"] = candidates.at[candidate,"length"] #TODO Add efficiency losses network.madd("Link", h2_links.index, bus0=h2_links.bus0.values + " H2", bus1=h2_links.bus1.values + " H2", p_min_pu=-1, p_nom_extendable=True, length=h2_links.length.values, capital_cost=costs.at['H2 pipeline','fixed']*h2_links.length.values, carrier="H2 pipeline", lifetime=costs.at['H2 pipeline','lifetime']) network.add("Carrier","battery") network.madd("Bus", nodes + " battery", location=nodes, carrier="battery") network.madd("Store", nodes + " battery", bus=nodes + " battery", e_cyclic=True, e_nom_extendable=True, carrier="battery", capital_cost=costs.at['battery storage','fixed'], lifetime=costs.at['battery storage','lifetime']) network.madd("Link", nodes + " battery charger", bus0=nodes, bus1=nodes + " battery", carrier="battery charger", efficiency=costs.at['battery inverter','efficiency']**0.5, capital_cost=costs.at['battery inverter','fixed'], p_nom_extendable=True, lifetime=costs.at['battery inverter','lifetime']) network.madd("Link", nodes + " battery discharger", bus0=nodes + " battery", bus1=nodes, carrier="battery discharger", efficiency=costs.at['battery inverter','efficiency']**0.5, marginal_cost=options['marginal_cost_storage'], p_nom_extendable=True, lifetime=costs.at['battery inverter','lifetime']) if options['methanation']: network.madd("Link", nodes + " Sabatier", bus0=nodes+" H2", bus1=["EU gas"]*len(nodes), bus2="co2 stored", p_nom_extendable=True, carrier="Sabatier", efficiency=costs.at["methanation","efficiency"], efficiency2=-costs.at["methanation","efficiency"]*costs.at['gas','CO2 intensity'], capital_cost=costs.at["methanation","fixed"], lifetime=costs.at['methanation','lifetime']) if options['helmeth']: network.madd("Link", nodes + " helmeth", bus0=nodes, bus1=["EU gas"]*len(nodes), bus2="co2 stored", carrier="helmeth", p_nom_extendable=True, efficiency=costs.at["helmeth","efficiency"], efficiency2=-costs.at["helmeth","efficiency"]*costs.at['gas','CO2 intensity'], capital_cost=costs.at["helmeth","fixed"], lifetime=costs.at['helmeth','lifetime']) if options['SMR']: network.madd("Link", nodes + " SMR CC", bus0=["EU gas"]*len(nodes), bus1=nodes+" H2", bus2="co2 atmosphere", bus3="co2 stored", p_nom_extendable=True, carrier="SMR CC", efficiency=costs.at["SMR CC","efficiency"], efficiency2=costs.at['gas','CO2 intensity']*(1-options["cc_fraction"]), efficiency3=costs.at['gas','CO2 intensity']*options["cc_fraction"], capital_cost=costs.at["SMR CC","fixed"], lifetime=costs.at['SMR CC','lifetime']) network.madd("Link", nodes + " SMR", bus0=["EU gas"]*len(nodes), bus1=nodes+" H2", bus2="co2 atmosphere", p_nom_extendable=True, carrier="SMR", efficiency=costs.at["SMR","efficiency"], efficiency2=costs.at['gas','CO2 intensity'], capital_cost=costs.at["SMR","fixed"], lifetime=costs.at['SMR','lifetime']) def add_land_transport(network): print("adding land transport") fuel_cell_share = get_parameter(options["land_transport_fuel_cell_share"]) electric_share = get_parameter(options["land_transport_electric_share"]) ice_share = 1 - fuel_cell_share - electric_share print("shares of FCEV, EV and ICEV are", fuel_cell_share, electric_share, ice_share) if ice_share < 0: print("Error, more FCEV and EV share than 1.") sys.exit() nodes = pop_layout.index if electric_share > 0: network.add("Carrier","Li ion") network.madd("Bus", nodes, location=nodes, suffix=" EV battery", carrier="Li ion") network.madd("Load", nodes, suffix=" land transport EV", bus=nodes + " EV battery", carrier="land transport EV", p_set=electric_share*(transport[nodes]+shift_df(transport[nodes],1)+shift_df(transport[nodes],2))/3.) p_nom = nodal_transport_data["number cars"]*0.011*electric_share #3-phase charger with 11 kW * x% of time grid-connected network.madd("Link", nodes, suffix= " BEV charger", bus0=nodes, bus1=nodes + " EV battery", p_nom=p_nom, carrier="BEV charger", p_max_pu=avail_profile[nodes], efficiency=0.9, #[B] #These were set non-zero to find LU infeasibility when availability = 0.25 #p_nom_extendable=True, #p_nom_min=p_nom, #capital_cost=1e6, #i.e. so high it only gets built where necessary ) if options["v2g"]: network.madd("Link", nodes, suffix=" V2G", bus1=nodes, bus0=nodes + " EV battery", p_nom=p_nom, carrier="V2G", p_max_pu=avail_profile[nodes], efficiency=0.9) #[B] if options["bev_dsm"]: network.madd("Store", nodes, suffix=" battery storage", bus=nodes + " EV battery", carrier="battery storage", e_cyclic=True, e_nom=nodal_transport_data["number cars"]*0.05*options["bev_availability"]*electric_share, #50 kWh battery http://www.zeit.de/mobilitaet/2014-10/auto-fahrzeug-bestand e_max_pu=1, e_min_pu=dsm_profile[nodes]) if fuel_cell_share > 0: network.madd("Load", nodes, suffix=" land transport fuel cell", bus=nodes + " H2", carrier="land transport fuel cell", p_set=fuel_cell_share/options['transport_fuel_cell_efficiency']*transport[nodes]) if ice_share > 0: network.madd("Load", nodes, suffix=" land transport oil", bus="EU oil", carrier="land transport oil", p_set=ice_share/options['transport_internal_combustion_efficiency']*transport[nodes]) def add_heat(network): print("adding heat") sectors = ["residential", "services"] nodes = create_nodes_for_heat_sector() #NB: must add costs of central heating afterwards (EUR 400 / kWpeak, 50a, 1% FOM from Fraunhofer ISE) urban_fraction = options['central_fraction']*pop_layout["urban"]/(pop_layout[["urban","rural"]].sum(axis=1)) # building retrofitting, exogenously reduce space heat demand if options["retrofitting"]["retro_exogen"]: dE = get_parameter(options["retrofitting"]["dE"]) print("retrofitting exogenously, assumed space heat reduction of ", dE) for sector in sectors: heat_demand[sector + " space"] = (1-dE)*heat_demand[sector + " space"] heat_systems = ["residential rural", "services rural", "residential urban decentral","services urban decentral", "urban central"] for name in heat_systems: name_type = "central" if name == "urban central" else "decentral" network.add("Carrier",name + " heat") network.madd("Bus", nodes[name] + " " + name + " heat", location=nodes[name], carrier=name + " heat") ## Add heat load for sector in sectors: if "rural" in name: factor = 1-urban_fraction[nodes[name]] elif "urban" in name: factor = urban_fraction[nodes[name]] else: factor = None if sector in name: heat_load = heat_demand[[sector + " water",sector + " space"]].groupby(level=1,axis=1).sum()[nodes[name]].multiply(factor) if name == "urban central": heat_load = heat_demand.groupby(level=1,axis=1).sum()[nodes[name]].multiply(urban_fraction[nodes[name]]*(1+options['district_heating_loss'])) network.madd("Load", nodes[name], suffix=" " + name + " heat", bus=nodes[name] + " " + name + " heat", carrier=name + " heat", p_set=heat_load) ## Add heat pumps heat_pump_type = "air" if "urban" in name else "ground" costs_name = "{} {}-sourced heat pump".format(name_type,heat_pump_type) cop = {"air" : ashp_cop, "ground" : gshp_cop} efficiency = cop[heat_pump_type][nodes[name]] if options["time_dep_hp_cop"] else costs.at[costs_name,'efficiency'] network.madd("Link", nodes[name], suffix=" {} {} heat pump".format(name,heat_pump_type), bus0=nodes[name], bus1=nodes[name] + " " + name + " heat", carrier="{} {} heat pump".format(name,heat_pump_type), efficiency=efficiency, capital_cost=costs.at[costs_name,'efficiency']*costs.at[costs_name,'fixed'], p_nom_extendable=True, lifetime=costs.at[costs_name,'lifetime']) if options["tes"]: network.add("Carrier",name + " water tanks") network.madd("Bus", nodes[name] + " " + name + " water tanks", location=nodes[name], carrier=name + " water tanks") network.madd("Link", nodes[name] + " " + name + " water tanks charger", bus0=nodes[name] + " " + name + " heat", bus1=nodes[name] + " " + name + " water tanks", efficiency=costs.at['water tank charger','efficiency'], carrier=name + " water tanks charger", p_nom_extendable=True) network.madd("Link", nodes[name] + " " + name + " water tanks discharger", bus0=nodes[name] + " " + name + " water tanks", bus1=nodes[name] + " " + name + " heat", carrier=name + " water tanks discharger", efficiency=costs.at['water tank discharger','efficiency'], p_nom_extendable=True) # [HP] 180 day time constant for centralised, 3 day for decentralised tes_time_constant_days = options["tes_tau"] if name_type == "decentral" else 180. network.madd("Store", nodes[name] + " " + name + " water tanks", bus=nodes[name] + " " + name + " water tanks", e_cyclic=True, e_nom_extendable=True, carrier=name + " water tanks", standing_loss=1-np.exp(-1/(24.*tes_time_constant_days)), capital_cost=costs.at[name_type + ' water tank storage','fixed']/(1.17e-3*40), #conversion from EUR/m^3 to EUR/MWh for 40 K diff and 1.17 kWh/m^3/K lifetime=costs.at[name_type + ' water tank storage','lifetime']) if options["boilers"]: network.madd("Link", nodes[name] + " " + name + " resistive heater", bus0=nodes[name], bus1=nodes[name] + " " + name + " heat", carrier=name + " resistive heater", efficiency=costs.at[name_type + ' resistive heater','efficiency'], capital_cost=costs.at[name_type + ' resistive heater','efficiency']*costs.at[name_type + ' resistive heater','fixed'], p_nom_extendable=True, lifetime=costs.at[name_type + ' resistive heater','lifetime']) network.madd("Link", nodes[name] + " " + name + " gas boiler", p_nom_extendable=True, bus0=["EU gas"]*len(nodes[name]), bus1=nodes[name] + " " + name + " heat", bus2="co2 atmosphere", carrier=name + " gas boiler", efficiency=costs.at[name_type + ' gas boiler','efficiency'], efficiency2=costs.at['gas','CO2 intensity'], capital_cost=costs.at[name_type + ' gas boiler','efficiency']*costs.at[name_type + ' gas boiler','fixed'], lifetime=costs.at[name_type + ' gas boiler','lifetime']) if options["solar_thermal"]: network.add("Carrier",name + " solar thermal") network.madd("Generator", nodes[name], suffix=" " + name + " solar thermal collector", bus=nodes[name] + " " + name + " heat", carrier=name + " solar thermal", p_nom_extendable=True, capital_cost=costs.at[name_type + ' solar thermal','fixed'], p_max_pu=solar_thermal[nodes[name]], lifetime=costs.at[name_type + ' solar thermal','lifetime']) if options["chp"]: if name == "urban central": #add gas CHP; biomass CHP is added in biomass section network.madd("Link", nodes[name] + " urban central gas CHP", bus0="EU gas", bus1=nodes[name], bus2=nodes[name] + " urban central heat", bus3="co2 atmosphere", carrier="urban central gas CHP", p_nom_extendable=True, capital_cost=costs.at['central gas CHP','fixed']*costs.at['central gas CHP','efficiency'], marginal_cost=costs.at['central gas CHP','VOM'], efficiency=costs.at['central gas CHP','efficiency'], efficiency2=costs.at['central gas CHP','efficiency']/costs.at['central gas CHP','c_b'], efficiency3=costs.at['gas','CO2 intensity'], lifetime=costs.at['central gas CHP','lifetime']) network.madd("Link", nodes[name] + " urban central gas CHP CC", bus0="EU gas", bus1=nodes[name], bus2=nodes[name] + " urban central heat", bus3="co2 atmosphere", bus4="co2 stored", carrier="urban central gas CHP CC", p_nom_extendable=True, capital_cost=costs.at['central gas CHP','fixed']*costs.at['central gas CHP','efficiency'] + costs.at['biomass CHP capture','fixed']*costs.at['gas','CO2 intensity'], marginal_cost=costs.at['central gas CHP','VOM'], efficiency=costs.at['central gas CHP','efficiency'] - costs.at['gas','CO2 intensity']*(costs.at['biomass CHP capture','electricity-input'] + costs.at['biomass CHP capture','compression-electricity-input']), efficiency2=costs.at['central gas CHP','efficiency']/costs.at['central gas CHP','c_b'] + costs.at['gas','CO2 intensity']*(costs.at['biomass CHP capture','heat-output'] + costs.at['biomass CHP capture','compression-heat-output'] - costs.at['biomass CHP capture','heat-output']), efficiency3=costs.at['gas','CO2 intensity']*(1-costs.at['biomass CHP capture','capture_rate']), efficiency4=costs.at['gas','CO2 intensity']*costs.at['biomass CHP capture','capture_rate'], lifetime=costs.at['central gas CHP','lifetime']) else: if options["micro_chp"]: network.madd("Link", nodes[name] + " " + name + " micro gas CHP", p_nom_extendable=True, bus0="EU gas", bus1=nodes[name], bus2=nodes[name] + " " + name + " heat", bus3="co2 atmosphere", carrier=name + " micro gas CHP", efficiency=costs.at['micro CHP','efficiency'], efficiency2=costs.at['micro CHP','efficiency-heat'], efficiency3=costs.at['gas','CO2 intensity'], capital_cost=costs.at['micro CHP','fixed'], lifetime=costs.at['micro CHP','lifetime']) if options['retrofitting']['retro_endogen']: print("adding retrofitting endogenously") # resample heat demand temporal 'heat_demand_r' depending on in config # specified temporal resolution, to not overestimate retrofitting hours = list(filter(re.compile(r'^\d+h$', re.IGNORECASE).search, opts)) if len(hours)==0: hours = [n.snapshots[1] - n.snapshots[0]] heat_demand_r = heat_demand.resample(hours[0]).mean() # retrofitting data 'retro_data' with 'costs' [EUR/m^2] and heat # demand 'dE' [per unit of original heat demand] for each country and # different retrofitting strengths [additional insulation thickness in m] retro_data = pd.read_csv(snakemake.input.retro_cost_energy, index_col=[0, 1], skipinitialspace=True, header=[0, 1]) # heated floor area [10^6 * m^2] per country floor_area = pd.read_csv(snakemake.input.floor_area, index_col=[0, 1]) network.add("Carrier", "retrofitting") # share of space heat demand 'w_space' of total heat demand w_space = {} for sector in sectors: w_space[sector] = heat_demand_r[sector + " space"] / \ (heat_demand_r[sector + " space"] + heat_demand_r[sector + " water"]) w_space["tot"] = ((heat_demand_r["services space"] + heat_demand_r["residential space"]) / heat_demand_r.groupby(level=[1], axis=1).sum()) for name in network.loads[network.loads.carrier.isin([x + " heat" for x in heat_systems])].index: node = network.buses.loc[name, "location"] ct = pop_layout.loc[node, "ct"] # weighting 'f' depending on the size of the population at the node f = urban_fraction[node] if "urban" in name else (1-urban_fraction[node]) if f == 0: continue # get sector name ("residential"/"services"/or both "tot" for urban central) sec = [x if x in name else "tot" for x in sectors][0] # get floor aread at node and region (urban/rural) in m^2 floor_area_node = ((pop_layout.loc[node].fraction * floor_area.loc[ct, "value"] * 10**6).loc[sec] * f) # total heat demand at node [MWh] demand = (network.loads_t.p_set[name].resample(hours[0]) .mean()) # space heat demand at node [MWh] space_heat_demand = demand * w_space[sec][node] # normed time profile of space heat demand 'space_pu' (values between 0-1), # p_max_pu/p_min_pu of retrofitting generators space_pu = (space_heat_demand / space_heat_demand.max()).to_frame(name=node) # minimum heat demand 'dE' after retrofitting in units of original heat demand (values between 0-1) dE = retro_data.loc[(ct, sec), ("dE")] # get addtional energy savings 'dE_diff' between the different retrofitting strengths/generators at one node dE_diff = abs(dE.diff()).fillna(1-dE.iloc[0]) # convert costs Euro/m^2 -> Euro/MWh capital_cost = retro_data.loc[(ct, sec), ("cost")] * floor_area_node / \ ((1 - dE) * space_heat_demand.max()) # number of possible retrofitting measures 'strengths' (set in list at config.yaml 'l_strength') # given in additional insulation thickness [m] # for each measure, a retrofitting generator is added at the node strengths = retro_data.columns.levels[1] # check that ambitious retrofitting has higher costs per MWh than moderate retrofitting if (capital_cost.diff() < 0).sum(): print( "warning, costs are not linear for ", ct, " ", sec) s = capital_cost[(capital_cost.diff() < 0)].index strengths = strengths.drop(s) # reindex normed time profile of space heat demand back to hourly resolution space_pu = (space_pu.reindex(index=heat_demand.index) .fillna(method="ffill")) # add for each retrofitting strength a generator with heat generation profile following the profile of the heat demand for strength in strengths: network.madd('Generator', [node], suffix=' retrofitting ' + strength + " " + name[6::], bus=name, carrier="retrofitting", p_nom_extendable=True, p_nom_max=dE_diff[strength] * space_heat_demand.max(), # maximum energy savings for this renovation strength p_max_pu=space_pu, p_min_pu=space_pu, country=ct, capital_cost=capital_cost[strength] * options['retrofitting']['cost_factor']) def create_nodes_for_heat_sector(): sectors = ["residential", "services"] # stores the different groups of nodes nodes = {} # rural are areas with low heating density and individual heating # urban are areas with high heating density # urban can be split into district heating (central) and individual heating (decentral) for sector in sectors: nodes[sector + " rural"] = pop_layout.index if options["central"]: urban_decentral_ct = pd.Index(["ES", "GR", "PT", "IT", "BG"]) nodes[sector + " urban decentral"] = pop_layout.index[pop_layout.ct.isin(urban_decentral_ct)] else: nodes[sector + " urban decentral"] = pop_layout.index # for central nodes, residential and services are aggregated nodes["urban central"] = pop_layout.index ^ nodes["residential urban decentral"] return nodes def add_biomass(network): print("adding biomass") nodes = pop_layout.index #biomass distributed at country level - i.e. transport within country allowed cts = pop_layout.ct.value_counts().index biomass_potentials = pd.read_csv(snakemake.input.biomass_potentials, index_col=0) network.add("Carrier","biogas") network.add("Carrier","solid biomass") network.madd("Bus", ["EU biogas"], location="EU", carrier="biogas") network.madd("Bus", ["EU solid biomass"], location="EU", carrier="solid biomass") network.madd("Store", ["EU biogas"], bus="EU biogas", carrier="biogas", e_nom=biomass_potentials.loc[cts,"biogas"].sum(), marginal_cost=costs.at['biogas','fuel'], e_initial=biomass_potentials.loc[cts,"biogas"].sum()) network.madd("Store", ["EU solid biomass"], bus="EU solid biomass", carrier="solid biomass", e_nom=biomass_potentials.loc[cts,"solid biomass"].sum(), marginal_cost=costs.at['solid biomass','fuel'], e_initial=biomass_potentials.loc[cts,"solid biomass"].sum()) network.madd("Link", ["biogas to gas"], bus0="EU biogas", bus1="EU gas", bus2="co2 atmosphere", carrier="biogas to gas", efficiency2=-costs.at['gas','CO2 intensity'], p_nom_extendable=True) #AC buses with district heating urban_central = network.buses.index[network.buses.carrier == "urban central heat"] if not urban_central.empty and options["chp"]: urban_central = urban_central.str[:-len(" urban central heat")] network.madd("Link", urban_central + " urban central solid biomass CHP", bus0="EU solid biomass", bus1=urban_central, bus2=urban_central + " urban central heat", carrier="urban central solid biomass CHP", p_nom_extendable=True, capital_cost=costs.at['central solid biomass CHP','fixed']*costs.at['central solid biomass CHP','efficiency'], marginal_cost=costs.at['central solid biomass CHP','VOM'], efficiency=costs.at['central solid biomass CHP','efficiency'], efficiency2=costs.at['central solid biomass CHP','efficiency-heat'], lifetime=costs.at['central solid biomass CHP','lifetime']) network.madd("Link", urban_central + " urban central solid biomass CHP CC", bus0="EU solid biomass", bus1=urban_central, bus2=urban_central + " urban central heat", bus3="co2 atmosphere", bus4="co2 stored", carrier="urban central solid biomass CHP CC", p_nom_extendable=True, capital_cost=costs.at['central solid biomass CHP','fixed']*costs.at['central solid biomass CHP','efficiency'] + costs.at['biomass CHP capture','fixed']*costs.at['solid biomass','CO2 intensity'], marginal_cost=costs.at['central solid biomass CHP','VOM'], efficiency=costs.at['central solid biomass CHP','efficiency'] - costs.at['solid biomass','CO2 intensity']*(costs.at['biomass CHP capture','electricity-input'] + costs.at['biomass CHP capture','compression-electricity-input']), efficiency2=costs.at['central solid biomass CHP','efficiency-heat'] + costs.at['solid biomass','CO2 intensity']*(costs.at['biomass CHP capture','heat-output'] + costs.at['biomass CHP capture','compression-heat-output'] - costs.at['biomass CHP capture','heat-output']), efficiency3=-costs.at['solid biomass','CO2 intensity']*costs.at['biomass CHP capture','capture_rate'], efficiency4=costs.at['solid biomass','CO2 intensity']*costs.at['biomass CHP capture','capture_rate'], lifetime=costs.at['central solid biomass CHP','lifetime']) def add_industry(network): print("adding industrial demand") nodes = pop_layout.index #1e6 to convert TWh to MWh industrial_demand = 1e6*pd.read_csv(snakemake.input.industrial_demand, index_col=0) solid_biomass_by_country = industrial_demand["solid biomass"].groupby(pop_layout.ct).sum() countries = solid_biomass_by_country.index network.madd("Bus", ["solid biomass for industry"], location="EU", carrier="solid biomass for industry") network.madd("Load", ["solid biomass for industry"], bus="solid biomass for industry", carrier="solid biomass for industry", p_set=solid_biomass_by_country.sum()/8760.) network.madd("Link", ["solid biomass for industry"], bus0="EU solid biomass", bus1="solid biomass for industry", carrier="solid biomass for industry", p_nom_extendable=True, efficiency=1.) network.madd("Link", ["solid biomass for industry CC"], bus0="EU solid biomass", bus1="solid biomass for industry", bus2="co2 atmosphere", bus3="co2 stored", carrier="solid biomass for industry CC", p_nom_extendable=True, capital_cost=costs.at["cement capture","fixed"]*costs.at['solid biomass','CO2 intensity'], efficiency=0.9, efficiency2=-costs.at['solid biomass','CO2 intensity']*costs.at["cement capture","capture_rate"], efficiency3=costs.at['solid biomass','CO2 intensity']*costs.at["cement capture","capture_rate"], lifetime=costs.at['cement capture','lifetime']) network.madd("Bus", ["gas for industry"], location="EU", carrier="gas for industry") network.madd("Load", ["gas for industry"], bus="gas for industry", carrier="gas for industry", p_set=industrial_demand.loc[nodes,"methane"].sum()/8760.) network.madd("Link", ["gas for industry"], bus0="EU gas", bus1="gas for industry", bus2="co2 atmosphere", carrier="gas for industry", p_nom_extendable=True, efficiency=1., efficiency2=costs.at['gas','CO2 intensity']) network.madd("Link", ["gas for industry CC"], bus0="EU gas", bus1="gas for industry", bus2="co2 atmosphere", bus3="co2 stored", carrier="gas for industry CC", p_nom_extendable=True, capital_cost=costs.at["cement capture","fixed"]*costs.at['gas','CO2 intensity'], efficiency=0.9, efficiency2=costs.at['gas','CO2 intensity']*(1-costs.at["cement capture","capture_rate"]), efficiency3=costs.at['gas','CO2 intensity']**costs.at["cement capture","capture_rate"], lifetime=costs.at['cement capture','lifetime']) network.madd("Load", nodes, suffix=" H2 for industry", bus=nodes + " H2", carrier="H2 for industry", p_set=industrial_demand.loc[nodes,"hydrogen"]/8760.) network.madd("Load", nodes, suffix=" H2 for shipping", bus=nodes + " H2", carrier="H2 for shipping", p_set = nodal_energy_totals.loc[nodes,["total international navigation","total domestic navigation"]].sum(axis=1)*1e6*options['shipping_average_efficiency']/costs.at["fuel cell","efficiency"]/8760.) if "EU oil" not in network.buses.index: network.madd("Bus", ["EU oil"], location="EU", carrier="oil") #use madd to get carrier inserted if "EU oil Store" not in network.stores.index: network.madd("Store", ["EU oil Store"], bus="EU oil", e_nom_extendable=True, e_cyclic=True, carrier="oil", capital_cost=0.) #could correct to e.g. 0.001 EUR/kWh * annuity and O&M if "EU oil" not in network.generators.index: network.add("Generator", "EU oil", bus="EU oil", p_nom_extendable=True, carrier="oil", capital_cost=0., marginal_cost=costs.at["oil",'fuel']) if options["oil_boilers"]: nodes_heat = create_nodes_for_heat_sector() for name in ["residential rural", "services rural", "residential urban decentral", "services urban decentral"]: network.madd("Link", nodes_heat[name] + " " + name + " oil boiler", p_nom_extendable=True, bus0="EU oil", bus1=nodes_heat[name] + " " + name + " heat", bus2="co2 atmosphere", carrier=name + " oil boiler", efficiency=costs.at['decentral oil boiler', 'efficiency'], efficiency2=costs.at['oil', 'CO2 intensity'], capital_cost=costs.at['decentral oil boiler', 'efficiency'] * costs.at[ 'decentral oil boiler', 'fixed'], lifetime=costs.at['decentral oil boiler','lifetime']) network.madd("Link", nodes + " Fischer-Tropsch", bus0=nodes + " H2", bus1="EU oil", bus2="co2 stored", carrier="Fischer-Tropsch", efficiency=costs.at["Fischer-Tropsch",'efficiency'], capital_cost=costs.at["Fischer-Tropsch",'fixed'], efficiency2=-costs.at["oil",'CO2 intensity']*costs.at["Fischer-Tropsch",'efficiency'], p_nom_extendable=True, lifetime=costs.at['Fischer-Tropsch','lifetime']) network.madd("Load", ["naphtha for industry"], bus="EU oil", carrier="naphtha for industry", p_set = industrial_demand.loc[nodes,"naphtha"].sum()/8760.) network.madd("Load", ["kerosene for aviation"], bus="EU oil", carrier="kerosene for aviation", p_set = nodal_energy_totals.loc[nodes,["total international aviation","total domestic aviation"]].sum(axis=1).sum()*1e6/8760.) #NB: CO2 gets released again to atmosphere when plastics decay or kerosene is burned #except for the process emissions when naphtha is used for petrochemicals, which can be captured with other industry process emissions #tco2 per hour co2 = network.loads.loc[["naphtha for industry","kerosene for aviation"],"p_set"].sum()*costs.at["oil",'CO2 intensity'] - industrial_demand.loc[nodes,"process emission from feedstock"].sum()/8760. network.madd("Load", ["oil emissions"], bus="co2 atmosphere", carrier="oil emissions", p_set=-co2) network.madd("Load", nodes, suffix=" low-temperature heat for industry", bus=[node + " urban central heat" if node + " urban central heat" in network.buses.index else node + " services urban decentral heat" for node in nodes], carrier="low-temperature heat for industry", p_set=industrial_demand.loc[nodes,"low-temperature heat"]/8760.) #remove today's industrial electricity demand by scaling down total electricity demand for ct in n.buses.country.unique(): loads = n.loads.index[(n.loads.index.str[:2] == ct) & (n.loads.carrier == "electricity")] factor = 1 - industrial_demand.loc[loads,"current electricity"].sum()/n.loads_t.p_set[loads].sum().sum() n.loads_t.p_set[loads] *= factor network.madd("Load", nodes, suffix=" industry electricity", bus=nodes, carrier="industry electricity", p_set=industrial_demand.loc[nodes,"electricity"]/8760.) network.madd("Bus", ["process emissions"], location="EU", carrier="process emissions") #this should be process emissions fossil+feedstock #then need load on atmosphere for feedstock emissions that are currently going to atmosphere via Link Fischer-Tropsch demand network.madd("Load", ["process emissions"], bus="process emissions", carrier="process emissions", p_set = -industrial_demand.loc[nodes,["process emission","process emission from feedstock"]].sum(axis=1).sum()/8760.) network.madd("Link", ["process emissions"], bus0="process emissions", bus1="co2 atmosphere", carrier="process emissions", p_nom_extendable=True, efficiency=1.) #assume enough local waste heat for CC network.madd("Link", ["process emissions CC"], bus0="process emissions", bus1="co2 atmosphere", bus2="co2 stored", carrier="process emissions CC", p_nom_extendable=True, capital_cost=costs.at["cement capture","fixed"], efficiency=(1-costs.at["cement capture","capture_rate"]), efficiency2=costs.at["cement capture","capture_rate"], lifetime=costs.at['cement capture','lifetime']) def add_waste_heat(network): print("adding possibility to use industrial waste heat in district heating") #AC buses with district heating urban_central = network.buses.index[network.buses.carrier == "urban central heat"] if not urban_central.empty: urban_central = urban_central.str[:-len(" urban central heat")] if options['use_fischer_tropsch_waste_heat']: network.links.loc[urban_central + " Fischer-Tropsch","bus3"] = urban_central + " urban central heat" network.links.loc[urban_central + " Fischer-Tropsch","efficiency3"] = 0.95 - network.links.loc[urban_central + " Fischer-Tropsch","efficiency"] if options['use_fuel_cell_waste_heat']: network.links.loc[urban_central + " H2 Fuel Cell","bus2"] = urban_central + " urban central heat" network.links.loc[urban_central + " H2 Fuel Cell","efficiency2"] = 0.95 - network.links.loc[urban_central + " H2 Fuel Cell","efficiency"] def decentral(n): n.lines.drop(n.lines.index,inplace=True) n.links.drop(n.links.index[n.links.carrier.isin(["DC","B2B"])],inplace=True) def remove_h2_network(n): nodes = pop_layout.index n.links.drop(n.links.index[n.links.carrier.isin(["H2 pipeline"])],inplace=True) n.stores.drop(["EU H2 Store"],inplace=True) if options['hydrogen_underground_storage']: h2_capital_cost = costs.at["gas storage","fixed"] #h2_capital_cost = costs.at["hydrogen underground storage","fixed"] else: h2_capital_cost = costs.at["hydrogen storage","fixed"] #put back nodal H2 storage n.madd("Store", nodes + " H2 Store", bus=nodes + " H2", e_nom_extendable=True, e_cyclic=True, carrier="H2 Store", capital_cost=h2_capital_cost) def get_parameter(item): """Check whether it depends on investment year""" if type(item) is dict: return item[investment_year] else: return item if __name__ == "__main__": # Detect running outside of snakemake and mock snakemake for testing if 'snakemake' not in globals(): from vresutils.snakemake import MockSnakemake snakemake = MockSnakemake( wildcards=dict(network='elec', simpl='', clusters='37', lv='1.0', opts='', planning_horizons='2020', sector_opts='120H-T-H-B-I-onwind+p3-dist1-cb48be3'), input=dict( network='../pypsa-eur/networks/{network}_s{simpl}_{clusters}_ec_lv{lv}_{opts}.nc', energy_totals_name='resources/energy_totals.csv', co2_totals_name='resources/co2_totals.csv', transport_name='resources/transport_data.csv', traffic_data = "data/emobility/", biomass_potentials='resources/biomass_potentials.csv', timezone_mappings='data/timezone_mappings.csv', heat_profile="data/heat_load_profile_BDEW.csv", costs="../technology-data/outputs/costs_{planning_horizons}.csv", h2_cavern = "data/hydrogen_salt_cavern_potentials.csv", profile_offwind_ac="../pypsa-eur/resources/profile_offwind-ac.nc", profile_offwind_dc="../pypsa-eur/resources/profile_offwind-dc.nc", busmap_s="../pypsa-eur/resources/busmap_{network}_s{simpl}.csv", busmap="../pypsa-eur/resources/busmap_{network}_s{simpl}_{clusters}.csv", clustered_pop_layout="resources/pop_layout_{network}_s{simpl}_{clusters}.csv", simplified_pop_layout="resources/pop_layout_{network}_s{simpl}.csv", industrial_demand="resources/industrial_energy_demand_{network}_s{simpl}_{clusters}.csv", heat_demand_urban="resources/heat_demand_urban_{network}_s{simpl}_{clusters}.nc", heat_demand_rural="resources/heat_demand_rural_{network}_s{simpl}_{clusters}.nc", heat_demand_total="resources/heat_demand_total_{network}_s{simpl}_{clusters}.nc", temp_soil_total="resources/temp_soil_total_{network}_s{simpl}_{clusters}.nc", temp_soil_rural="resources/temp_soil_rural_{network}_s{simpl}_{clusters}.nc", temp_soil_urban="resources/temp_soil_urban_{network}_s{simpl}_{clusters}.nc", temp_air_total="resources/temp_air_total_{network}_s{simpl}_{clusters}.nc", temp_air_rural="resources/temp_air_rural_{network}_s{simpl}_{clusters}.nc", temp_air_urban="resources/temp_air_urban_{network}_s{simpl}_{clusters}.nc", cop_soil_total="resources/cop_soil_total_{network}_s{simpl}_{clusters}.nc", cop_soil_rural="resources/cop_soil_rural_{network}_s{simpl}_{clusters}.nc", cop_soil_urban="resources/cop_soil_urban_{network}_s{simpl}_{clusters}.nc", cop_air_total="resources/cop_air_total_{network}_s{simpl}_{clusters}.nc", cop_air_rural="resources/cop_air_rural_{network}_s{simpl}_{clusters}.nc", cop_air_urban="resources/cop_air_urban_{network}_s{simpl}_{clusters}.nc", solar_thermal_total="resources/solar_thermal_total_{network}_s{simpl}_{clusters}.nc", solar_thermal_urban="resources/solar_thermal_urban_{network}_s{simpl}_{clusters}.nc", solar_thermal_rural="resources/solar_thermal_rural_{network}_s{simpl}_{clusters}.nc", retro_cost_energy = "resources/retro_cost_{network}_s{simpl}_{clusters}.csv", floor_area = "resources/floor_area_{network}_s{simpl}_{clusters}.csv" ), output=['results/version-cb48be3/prenetworks/{network}_s{simpl}_{clusters}_lv{lv}__{sector_opts}_{planning_horizons}.nc'] ) import yaml with open('config.yaml', encoding='utf8') as f: snakemake.config = yaml.safe_load(f) logging.basicConfig(level=snakemake.config['logging_level']) timezone_mappings = pd.read_csv(snakemake.input.timezone_mappings,index_col=0,squeeze=True,header=None) options = snakemake.config["sector"] opts = snakemake.wildcards.sector_opts.split('-') investment_year=int(snakemake.wildcards.planning_horizons[-4:]) n = pypsa.Network(snakemake.input.network, override_component_attrs=override_component_attrs) Nyears = n.snapshot_weightings.sum()/8760. pop_layout = pd.read_csv(snakemake.input.clustered_pop_layout,index_col=0) pop_layout["ct"] = pop_layout.index.str[:2] ct_total = pop_layout.total.groupby(pop_layout["ct"]).sum() pop_layout["ct_total"] = pop_layout["ct"].map(ct_total.get) pop_layout["fraction"] = pop_layout["total"]/pop_layout["ct_total"] simplified_pop_layout = pd.read_csv(snakemake.input.simplified_pop_layout,index_col=0) costs = prepare_costs(snakemake.input.costs, snakemake.config['costs']['USD2013_to_EUR2013'], snakemake.config['costs']['discountrate'], Nyears, snakemake.config['costs']['lifetime']) remove_elec_base_techs(n) n.loads["carrier"] = "electricity" n.buses["location"] = n.buses.index update_wind_solar_costs(n, costs) if snakemake.config["foresight"]=='myopic': add_lifetime_wind_solar(n) add_carrier_buses(n,snakemake.config['existing_capacities']['conventional_carriers']) add_co2_tracking(n) add_generation(n) add_storage(n) for o in opts: if o[:4] == "wave": wave_cost_factor = float(o[4:].replace("p",".").replace("m","-")) print("Including wave generators with cost factor of", wave_cost_factor) add_wave(n, wave_cost_factor) if o[:4] == "dist": snakemake.config["sector"]['electricity_distribution_grid'] = True snakemake.config["sector"]['electricity_distribution_grid_cost_factor'] = float(o[4:].replace("p",".").replace("m","-")) nodal_energy_totals, heat_demand, ashp_cop, gshp_cop, solar_thermal, transport, avail_profile, dsm_profile, co2_totals, nodal_transport_data = prepare_data(n) if "nodistrict" in opts: options["central"] = False if "T" in opts: add_land_transport(n) if "H" in opts: add_heat(n) if "B" in opts: add_biomass(n) if "I" in opts: add_industry(n) if "I" in opts and "H" in opts: add_waste_heat(n) if options['dac']: add_dac(n) if "decentral" in opts: decentral(n) if "noH2network" in opts: remove_h2_network(n) for o in opts: m = re.match(r'^\d+h$', o, re.IGNORECASE) if m is not None: n = average_every_nhours(n, m.group(0)) break else: logger.info("No resampling") #process CO2 limit limit = get_parameter(snakemake.config["co2_budget"]) print("CO2 limit set to",limit) for o in opts: if "cb" in o: path_cb = snakemake.config['results_dir'] + snakemake.config['run'] + '/csvs/' if not os.path.exists(path_cb): os.makedirs(path_cb) try: CO2_CAP=pd.read_csv(path_cb + 'carbon_budget_distribution.csv', index_col=0) except: build_carbon_budget(o) CO2_CAP=pd.read_csv(path_cb + 'carbon_budget_distribution.csv', index_col=0) limit=CO2_CAP.loc[investment_year] print("overriding CO2 limit with scenario limit",limit) for o in opts: if "Co2L" in o: limit = o[o.find("Co2L")+4:] limit = float(limit.replace("p",".").replace("m","-")) print("overriding CO2 limit with scenario limit",limit) print("adding CO2 budget limit as per unit of 1990 levels of",limit) add_co2limit(n, Nyears, limit) for o in opts: if o[:10] == 'linemaxext': maxext = float(o[10:])*1e3 print("limiting new HVAC and HVDC extensions to",maxext,"MW") n.lines['s_nom_max'] = n.lines['s_nom'] + maxext hvdc = n.links.index[n.links.carrier == 'DC'] n.links.loc[hvdc,'p_nom_max'] = n.links.loc[hvdc,'p_nom'] + maxext if snakemake.config["sector"]['electricity_distribution_grid']: insert_electricity_distribution_grid(n) for o in opts: if "+" in o: oo = o.split("+") carrier_list=np.hstack((n.generators.carrier.unique(), n.links.carrier.unique(), n.stores.carrier.unique(), n.storage_units.carrier.unique())) suptechs = map(lambda c: c.split("-", 2)[0], carrier_list) if oo[0].startswith(tuple(suptechs)): carrier = oo[0] attr_lookup = {"p": "p_nom_max", "c": "capital_cost"} attr = attr_lookup[oo[1][0]] factor = float(oo[1][1:]) #beware if factor is 0 and p_nom_max is np.inf, 0*np.inf is nan if carrier == "AC": # lines do not have carrier n.lines[attr] *= factor else: comps = {"Generator", "Link", "StorageUnit"} if attr=='p_nom_max' else {"Generator", "Link", "StorageUnit", "Store"} for c in n.iterate_components(comps): if carrier=='solar': sel = c.df.carrier.str.contains(carrier) & ~c.df.carrier.str.contains("solar rooftop") else: sel = c.df.carrier.str.contains(carrier) c.df.loc[sel,attr] *= factor print("changing", attr ,"for",carrier,"by factor",factor) if snakemake.config["sector"]['gas_distribution_grid']: insert_gas_distribution_costs(n) if snakemake.config["sector"]['electricity_grid_connection']: add_electricity_grid_connection(n) n.export_to_netcdf(snakemake.output[0])