Merge pull request #74 from PyPSA/retro

Retro

Snakefile

@@ -278,6 +278,23 @@ rule build_industrial_demand:
     resources: mem_mb=1000
     script: 'scripts/build_industrial_demand.py'
 
+if config['sector']['retrofitting'].get('retro_endogen', True):
+	rule build_retro_cost:
+	    input:
+    		building_stock="data/retro/data_building_stock.csv",
+    		u_values_PL="data/retro/u_values_poland.csv",
+    		tax_w="data/retro/electricity_taxes_eu.csv",
+    		construction_index="data/retro/comparative_level_investment.csv",
+    		average_surface="data/retro/average_surface_components.csv",
+    		floor_area_missing="data/retro/floor_area_missing.csv",
+    		clustered_pop_layout="resources/pop_layout_{network}_s{simpl}_{clusters}.csv",
+    		cost_germany="data/retro/retro_cost_germany.csv",
+    		window_assumptions="data/retro/window_assumptions.csv"
+	    output:
+    		retro_cost="resources/retro_cost_{network}_s{simpl}_{clusters}.csv",
+    		floor_area="resources/floor_area_{network}_s{simpl}_{clusters}.csv"
+	    script: "scripts/build_retro_cost.py"
+
 
 rule prepare_sector_network:
     input:
@@ -285,6 +302,7 @@ rule prepare_sector_network:
         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",
@@ -314,7 +332,9 @@ rule prepare_sector_network:
         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"
+        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: config['results_dir']  +  config['run'] + '/prenetworks/{network}_s{simpl}_{clusters}_lv{lv}_{opts}_{sector_opts}_{planning_horizons}.nc'
     threads: 1
     resources: mem_mb=2000

config.default.yaml

@@ -101,10 +101,16 @@ sector:
   'shipping_average_efficiency' : 0.4 #For conversion of fuel oil to propulsion in 2011
   'time_dep_hp_cop' : True
   'space_heating_fraction' : 1.0 #fraction of space heating active
-  'retrofitting' : False
+  'retrofitting' :
-  'retroI-fraction' : 0.25
+      'retro_exogen': False # space heat demand savings exogenously
-  'retroII-fraction' : 0.55
+      'dE': 0.4            # reduction of space heat demand (applied before losses in DH)
-  'retrofitting-cost_factor' : 1.0
+      'retro_endogen': True  # co-optimise space heat savings
+      'cost_factor' : 1.0
+      'interest_rate': 0.04  # for investment in building components
+      'annualise_cost': True  # annualise the investment costs
+      'tax_weighting': False   # weight costs depending on taxes in countries
+      'construction_index': True   # weight costs depending on labour/material costs per ct
+      'l_strength': ["0.076", "0.197"] # additional insulation thickness[m], determines number of retro steps(=generators per bus) and maximum possible savings
   'tes' : True
   'tes_tau' : 3.
   'boilers' : True
@@ -130,7 +136,6 @@ sector:
   'gas_distribution_grid_cost_factor' : 1.0  #multiplies cost in data/costs.csv
 
 costs:
-  year: 2030
   lifetime: 25 #default lifetime
   # From a Lion Hirth paper, also reflects average of Noothout et al 2016
   discountrate: 0.07

data/retro/average_surface_components.csv

@@ -0,0 +1,7 @@
+,Dwelling,Ceilling,Standard component surfaces (m2),component,surfaces,(m2),,
+Building type,Space(m²),Height(m),Roof,Facade,Floor,Windows,,
+Single/two family house,120,2.5,90,166,63,29,,
+Large apartment house,1457,2.5,354,1189,354,380,,
+Apartment house,5276,,598.337,2992.1,598.337,756,tabula ,http://webtool.building-typology.eu/#pdfes
+,,,,,,,,
+"Source: https://link.springer.com/article/10.1007/s12053-010-9090-6  ,p.4",,,,,,,,

data/retro/comparative_level_investment.csv

@@ -0,0 +1,49 @@
+NA_ITEM,Price level indices (EU28=100),,,,,,,,,
+PPP_CAT,Actual individual consumption,,,,,,,,,
+,,,,,,,,,,
+GEO/TIME,2009,2010,2011,2012,2013,2014,2015,2016,2017,2018
+European Union - 28 countries,100.0,100.0,100.0,100.0,100.0,100.0,100.0,100.0,100.0,100.0
+Belgium,113.6,111.9,112.4,111.5,111.0,108.9,106.3,110.3,112.3,112.5
+Bulgaria,47.1,45.7,45.5,45.0,44.2,42.6,42.2,43.2,45.1,46.3
+Czech Republic,64.5,66.6,68.9,66.9,63.3,58.3,58.4,60.5,62.4,65.0
+Denmark,141.7,140.0,139.9,140.0,139.3,138.5,135.0,140.0,138.9,138.1
+Germany,104.6,103.1,102.2,101.1,102.5,101.5,100.4,102.6,103.7,104.1
+Estonia,67.5,66.0,67.2,67.6,69.9,69.9,68.9,71.0,73.9,76.3
+Ireland,129.9,122.7,122.5,120.5,123.2,124.9,122.2,126.5,129.1,129.2
+Greece,93.6,95.4,94.9,91.9,87.8,83.8,81.0,82.3,83.0,81.8
+Spain,97.5,98.7,98.5,95.8,95.1,92.7,90.0,92.7,93.7,93.7
+France,111.2,109.9,109.6,108.7,107.0,106.0,104.0,105.8,107.1,107.4
+Croatia,70.2,70.1,68.1,65.5,64.5,62.5,60.7,61.3,63.0,64.0
+Italy,103.6,100.4,101.5,101.1,102.3,102.6,100.3,101.1,101.6,101.4
+Cyprus,92.0,94.6,95.8,96.0,95.2,92.0,88.5,89.8,91.2,90.6
+Latvia,68.1,62.3,65.5,65.9,66.0,66.0,64.2,66.9,68.3,69.5
+Lithuania,60.3,57.8,58.3,58.0,57.8,56.9,55.9,58.3,60.0,61.4
+Luxembourg,130.0,136.5,136.0,135.8,135.1,135.7,132.1,137.0,139.9,141.6
+Hungary,58.2,57.4,56.4,54.9,54.4,53.4,53.3,56.2,59.4,59.0
+Malta,75.8,76.6,78.0,78.0,80.8,80.5,79.8,81.4,81.9,83.4
+Netherlands,108.5,112.3,112.7,111.3,111.9,111.9,109.6,113.8,114.6,114.8
+Austria,109.9,109.2,110.1,108.9,109.1,109.1,107.2,110.2,112.8,113.7
+Poland,53.1,55.2,53.7,52.1,52.4,52.5,51.1,50.9,53.5,54.3
+Portugal,85.2,85.0,85.3,82.7,81.1,80.4,78.7,81.6,83.5,84.6
+Romania,49.1,46.9,47.7,45.6,47.8,47.6,47.2,46.8,48.0,48.6
+Slovenia,85.3,84.3,83.7,81.8,82.1,81.5,79.8,82.3,82.7,83.8
+Slovakia,66.6,62.5,63.4,63.4,63.4,63.3,62.3,63.6,65.4,66.1
+Finland,121.0,120.3,121.6,121.8,124.0,122.9,119.6,122.8,123.3,123.4
+Sweden,109.5,124.6,131.7,134.3,140.5,133.6,128.8,135.3,134.5,126.9
+United Kingdom,107.5,111.4,111.3,118.6,117.0,123.6,134.7,123.5,117.6,117.7
+Iceland,94.9,107.6,109.6,111.6,116.0,123.4,132.5,154.5,172.3,163.7
+Norway,142.4,158.8,165.3,172.5,166.9,157.2,152.2,155.0,157.3,155.4
+Switzerland,131.6,146.4,161.7,160.6,155.1,153.0,167.0,169.8,167.1,159.1
+Candidate and potential candidate countries except Turkey and Kosovo (under United Nations Security Council Resolution 1244/99),48.0,45.6,47.1,44.8,46.4,45.2,43.4,44.4,46.0,47.5
+Montenegro,52.3,49.5,49.3,50.1,50.5,49.3,48.0,48.7,50.5,51.1
+North Macedonia,41.4,41.3,42.7,42.1,42.5,41.9,40.9,41.7,43.2,43.3
+Albania,46.2,42.8,42.1,40.6,41.9,41.5,39.8,43.0,43.5,46.6
+Serbia,48.3,45.0,48.0,44.5,47.3,45.5,43.1,43.8,46.1,47.9
+Turkey,55.4,61.2,54.7,58.5,57.7,51.6,50.5,50.2,45.4,37.0
+Bosnia and Herzegovina,51.6,50.7,50.6,49.2,49.1,48.4,47.0,47.5,48.2,48.9
+Kosovo (under United Nations Security Council Resolution 1244/99),:,:,:,:,:,:,:,:,:,:
+United States,92.4,98,93.3,101.2,100.3,99,115.9,121.1,120.8,115.2
+Japan,115.1,126.1,127.8,133.8,101.7,94.8,96.5,113,109.4,103.9
+,,,,,,,,,,
+"Source: Eurostat Purchasing power parities (PPPs), price level indices and real expenditures for ESA 2010 aggregates (2019)",,,,,,,,,,
+https://ec.europa.eu/eurostat/statistics-explained/index.php?title=Comparative_price_levels_for_investment,,,,,,,,,,

data/retro/electricity_taxes_eu.csv

@@ -0,0 +1,164 @@
+Electricity prices for household consumers - bi-annual data (from 2007 onwards) [nrg_pc_204],,,,
+,,,,
+Last update,30.10.19,,,
+Extracted on,14.11.19,,,
+Source of data,Eurostat,,,
+,,,,
+PRODUCT,Electrical energy,,,
+CONSOM,Band DC : 2 500 kWh < Consumption < 5 000 kWh,,,
+UNIT,Kilowatt-hour,,,
+TIME,2018S1,,,
+,,,,
+CURRENCY,Euro,Euro,Euro,
+GEO/TAX,Excluding taxes and levies,Excluding VAT and other recoverable taxes and levies,All taxes and levies included,% cost without taxes
+European Union - 28 countries,0.1285,0.1756,0.2052,0.626218323586745
+"Euro area (EA11-2000, EA12-2006, EA13-2007, EA15-2008, EA16-2010, EA17-2013, EA18-2014, EA19)",0.1331,0.1855,0.2188,0.608318098720293
+Belgium,0.1903,0.2279,0.2733,0.696304427369191
+Bulgaria,0.0816,0.0816,0.0979,0.833503575076609
+Czech Republic,0.1286,0.1298,0.1573,0.817546090273363
+Denmark,0.1011,0.2501,0.3126,0.32341650671785
+Germany,0.1379,0.2510,0.2987,0.461667224640107
+Estonia,0.0989,0.1123,0.1348,0.733679525222552
+Ireland,0.1846,0.2087,0.2369,0.779231743351625
+Greece,0.1132,0.1482,0.1672,0.677033492822967
+Spain,0.1873,0.1969,0.2383,0.785984053713806
+France,0.1134,0.1492,0.1748,0.648741418764302
+Croatia,0.1020,0.1160,0.1311,0.778032036613272
+Italy,0.1285,0.1873,0.2067,0.621673923560716
+Cyprus,0.1445,0.1606,0.1893,0.763338615953513
+Latvia,0.1035,0.1266,0.1531,0.676028739386022
+Lithuania,0.0771,0.0906,0.1097,0.702825888787603
+Luxembourg,0.1283,0.1547,0.1671,0.767803710353082
+Hungary,0.0885,0.0885,0.1123,0.78806767586821
+Malta,0.1209,0.1224,0.1285,0.940856031128405
+Netherlands,0.1187,0.1410,0.1706,0.6957796014068
+Austria,0.1232,0.1638,0.1966,0.626653102746694
+Poland,0.0906,0.1146,0.1410,0.642553191489362
+Portugal,0.1007,0.1826,0.2246,0.448352626892253
+Romania,0.0990,0.1120,0.1333,0.742685671417854
+Slovenia,0.1108,0.1322,0.1613,0.686918784872908
+Slovakia,0.0942,0.1305,0.1566,0.601532567049808
+Finland,0.1074,0.1300,0.1612,0.666253101736973
+Sweden,0.1202,0.1513,0.1891,0.635642517186674
+United Kingdom,0.1347,0.1797,0.1887,0.713831478537361
+Iceland,0.1222,0.1246,0.1545,0.790938511326861
+Liechtenstein,:,:,:,#VALUE!
+Norway,0.1254,0.1434,0.1751,0.716162193032553
+Montenegro,0.0828,0.0844,0.1024,0.80859375
+North Macedonia,0.0662,0.0662,0.0781,0.847631241997439
+Albania,:,:,:,#VALUE!
+Serbia,0.0539,0.0587,0.0705,0.764539007092199
+Turkey,0.0727,0.0766,0.0904,0.804203539823009
+Bosnia and Herzegovina,0.0722,0.0738,0.0864,0.835648148148148
+Kosovo (under United Nations Security Council Resolution 1244/99),0.0569,0.0586,0.0633,0.898894154818325
+Moldova,0.1020,0.1020,0.1020,1
+Ukraine,0.0342,0.0342,0.0410,0.834146341463415
+,,,0.157271052631579,
+Special value:,,,,
+:,not available,,,
+,,,,
+PRODUCT,Electrical energy,,,
+CONSOM,Band DC : 2 500 kWh < Consumption < 5 000 kWh,,,
+UNIT,Kilowatt-hour,,,
+TIME,2018S2,,,
+,,,,
+CURRENCY,Euro,Euro,Euro,
+GEO/TAX,Excluding taxes and levies,Excluding VAT and other recoverable taxes and levies,All taxes and levies included,
+European Union - 28 countries,0.1329,0.1810,0.2113,
+"Euro area (EA11-2000, EA12-2006, EA13-2007, EA15-2008, EA16-2010, EA17-2013, EA18-2014, EA19)",0.1376,0.1902,0.2242,
+Belgium,0.1998,0.2429,0.2937,
+Bulgaria,0.0838,0.0838,0.1005,
+Czechia,0.1299,0.1311,0.1586,
+Denmark,0.1116,0.2499,0.3123,
+Germany (until 1990 former territory of the FRG),0.1378,0.2521,0.3000,
+Estonia,0.1048,0.1182,0.1418,
+Ireland,0.2006,0.2237,0.2539,
+Greece,0.1125,0.1458,0.1646,
+Spain,0.1947,0.2047,0.2477,
+France,0.1168,0.1537,0.1799,
+Croatia,0.1028,0.1169,0.1321,
+Italy,0.1416,0.1964,0.2161,
+Cyprus,0.1745,0.1850,0.2183,
+Latvia,0.1041,0.1249,0.1511,
+Lithuania,0.0771,0.0906,0.1097,
+Luxembourg,0.1302,0.1566,0.1691,
+Hungary,0.0880,0.0880,0.1118,
+Malta,0.1229,0.1244,0.1306,
+Netherlands,0.1212,0.1420,0.1707,
+Austria,0.1265,0.1676,0.2012,
+Poland,0.0889,0.1135,0.1396,
+Portugal,0.1028,0.1864,0.2293,
+Romania,0.0964,0.1107,0.1317,
+Slovenia,0.1125,0.1342,0.1638,
+Slovakia,0.0849,0.1218,0.1462,
+Finland,0.1144,0.1369,0.1698,
+Sweden,0.1287,0.1592,0.1990,
+United Kingdom,0.1401,0.1927,0.2024,
+Iceland,0.1152,0.1175,0.1457,
+Liechtenstein,:,:,:,
+Norway,0.1382,0.1562,0.1907,
+Montenegro,0.0829,0.0848,0.1030,
+North Macedonia,0.0667,0.0667,0.0787,
+Albania,0.0759,0.0759,0.0910,
+Serbia,0.0542,0.0591,0.0709,
+Turkey,0.0688,0.0726,0.0857,
+Bosnia and Herzegovina,0.0729,0.0744,0.0871,
+Kosovo (under United Nations Security Council Resolution 1244/99),0.0579,0.0591,0.0638,
+Moldova,0.0960,0.0960,0.1029,
+Ukraine,0.0342,0.0342,0.0410,
+,,,,
+Special value:,,,,
+:,not available,,,
+,,,,
+PRODUCT,Electrical energy,,,
+CONSOM,Band DC : 2 500 kWh < Consumption < 5 000 kWh,,,
+UNIT,Kilowatt-hour,,,
+TIME,2019S1,,,
+,,,,
+CURRENCY,Euro,Euro,Euro,
+GEO/TAX,Excluding taxes and levies,Excluding VAT and other recoverable taxes and levies,All taxes and levies included,
+European Union - 28 countries,0.1351,0.1841,0.2147,
+"Euro area (EA11-2000, EA12-2006, EA13-2007, EA15-2008, EA16-2010, EA17-2013, EA18-2014, EA19)",0.1396,0.1928,0.2270,
+Belgium,0.1965,0.2355,0.2839,
+Bulgaria,0.0831,0.0831,0.0997,
+Czechia,0.1433,0.1444,0.1748,
+Denmark,0.1084,0.2387,0.2984,
+Germany (until 1990 former territory of the FRG),0.1473,0.2595,0.3088,
+Estonia,0.0982,0.1131,0.1357,
+Ireland,0.2027,0.2134,0.2423,
+Greece,0.1139,0.1482,0.1650,
+Spain,0.1889,0.1986,0.2403,
+France,0.1138,0.1508,0.1765,
+Croatia,0.1028,0.1169,0.1321,
+Italy,0.1432,0.2090,0.2301,
+Cyprus,0.1762,0.1867,0.2203,
+Latvia,0.1136,0.1347,0.1629,
+Lithuania,0.0947,0.1037,0.1255,
+Luxembourg,0.1326,0.1666,0.1798,
+Hungary,0.0882,0.0882,0.1120,
+Malta,0.1228,0.1243,0.1305,
+Netherlands,0.1357,0.1708,0.2052,
+Austria,0.1316,0.1695,0.2034,
+Poland,0.0884,0.1092,0.1343,
+Portugal,0.1103,0.1751,0.2154,
+Romania,0.0983,0.1141,0.1358,
+Slovenia,0.1125,0.1339,0.1634,
+Slovakia,0.0962,0.1314,0.1577,
+Finland,0.1173,0.1398,0.1734,
+Sweden,0.1297,0.1612,0.2015,
+United Kingdom,0.1450,0.2021,0.2122,
+Iceland,0.1112,0.1134,0.1406,
+Liechtenstein,:,:,:,
+Norway,0.1360,0.1529,0.1867,
+Montenegro,0.0834,0.0850,0.1032,
+North Macedonia,:,:,:,
+Albania,:,:,:,
+Serbia,0.0541,0.0589,0.0706,
+Turkey,0.0684,0.0718,0.0847,
+Bosnia and Herzegovina,0.0729,0.0746,0.0873,
+Kosovo (under United Nations Security Council Resolution 1244/99),0.0537,0.0556,0.0600,
+Moldova,0.0936,0.0936,0.0936,
+Ukraine,0.0369,0.0369,0.0442,
+,,,,
+Special value:,,,,
+:,not available,,,

data/retro/floor_area_missing.csv

@@ -0,0 +1,17 @@
+country,sector,estimated,value,source,,comments,population [in Million],
+AL,residential,0,64,p.13 1.6 million m² = 2.5% of total floor area,https://www.buildup.eu/sites/default/files/content/sled_albania_residential_building_eng.pdf,,,
+AL,services,0,,,,,,
+BA,residential,0,125.89,Tabula,https://episcope.eu/building-typology/country/ba/,strong differences ? other source claims more than 300 Million m²,,https://www.buildup.eu/sites/default/files/content/sled_serbia_building_eng.pdf
+BA,services,0,,,,,,
+RS,residential,0,72.3,Odyssee(2011),https://odyssee.enerdata.net/database/,,,
+RS,services,0,,,,,,
+MK,residential,0,,"Worldbank p.7 Skopje 75% residential, 25% commercial",http://documents.albankaldawli.org/curated/ar/838951574180734318/pdf/Project-Information-Document-North-Macedonia-Public-Sector-Energy-Efficiency-Project-P149990.pdf,15 % live in illegal constructed buildings ? not part of the statistics,2.1,
+MK,services,0,,,,,,
+ME,residential,0,19.625,p.13 0.314 million m² = 1.6% of total floor area,buildup.eu/sites/default/files/content/sled_montenegro_building_eng.pdf,Only 50 % of the floor area is heated p.12,,buildup.eu/sites/default/files/content/sled_montenegro_building_eng.pdf
+ME,services,0,,,,,,
+CH,residential,0,99.45,Odyssee(2015),,,,
+CH,services,1,78.1392857142857,p.8 44%floor area is services,https://bta.climate-kic.org/wp-content/uploads/2018/04/171123-CK-BTA-DEF-BMB_SWITZERLAND_.pdf,,,
+NO,residential,0,121.55,Odyssee(2015),,,,
+NO,services,0,115.21,Odyssee(2015),,,,
+PL,residential,0,1028.41,EU Building Database,,,,
+PL,services,0,498.84,EU Building Database,,,,

data/retro/retro_cost_germany.csv

@@ -0,0 +1,7 @@
+component,cost_fix,cost_var,life_time,comment,additional source
+wall,70.34,2.36,40,Agora Energiewende p.110,
+floor,39.39,1.3,40,Agora Energiewende p.110,
+roof,75.61,1.3,40,Agora Energiewende p.110,https://www.baulinks.de/webplugin/2018/1524.php4
+window,nan,nan,35,,
+source: p.37 https://www.umweltbundesamt.de/sites/default/files/medien/1410/publikationen/2019-10-29_texte_132-2019_energieaufwand-gebaeudekonzepte.pdf,,,https://www.agora-energiewende.de/en/publications/building-sector-efficiency-a-crucial-component-of-the-energy-transition/,,
+,,,p.115,,

data/retro/u_values_poland.csv

@@ -0,0 +1,9 @@
+component,Before 1945,1945 - 1969,1970 - 1979,1980 - 1989,1990 - 1999,2000 - 2010,Post 2010,sector
+Walls,1.7,1.4,0.9,0.9,0.6,0.4,1.7,residential
+Windows,4.6,3.6,2.6,2.6,2.1,2.1,2.1,residential
+Roof,0.8,0.7,0.6,0.6,0.6,0.4,0.33,residential
+Floor,1.9,1.4,1.2,1.1,0.9,0.6,0.45,residential
+Walls,1.3,1.3,1.3,0.8,0.6,0.6,0.6,services
+Windows,4.7,3.7,2.6,2.6,2.3,2.1,2.1,services
+Roof,1,0.9,0.7,0.5,0.3,0.3,0.3,services
+Floor,1.6,1.2,1.2,1.1,1,0.7,0.7,services

data/retro/window_assumptions.csv

@@ -0,0 +1,8 @@
+strength,u_value,cost,u_limit,comment
+[m],[W/m^2K],EUR/m^2,[W/m^2K],
+0.076,1.34,180.08,3.5,Double-glazing
+0.197,0.8,225,1.3,Triple-glazing
+,,,,
+"source: https://www.agora-energiewende.de/en/publications/building-sector-efficiency-a-crucial-component-of-the-energy-transition/
+p.115
+",,,,

scripts/add_existing_baseyear.py

@@ -442,7 +442,8 @@ if __name__ == "__main__":
     costs = prepare_costs(snakemake.input.costs,
                           snakemake.config['costs']['USD2013_to_EUR2013'],
                           snakemake.config['costs']['discountrate'],
-                          Nyears)
+                          Nyears,
+                          snakemake.config['costs']['lifetime'])
 
     grouping_years=snakemake.config['existing_capacities']['grouping_years']
     add_power_capacities_installed_before_baseyear(n, grouping_years, costs, baseyear)

scripts/build_retro_cost.py

@@ -0,0 +1,477 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Jan 20 14:57:21 2020
+
+@author: bw0928
+
+*****************************************************************************
+This script calculates cost-energy_saving-curves for retrofitting
+for the EU-37 countries, based on the building stock data from hotmaps and
+the EU building stock database
+*****************************************************************************
+
+Structure:
+
+    (1) set assumptions and parameters
+    (2) read and prepare data
+    (3) calculate (€-dE-curves)
+    (4) save in csv
+
+*****************************************************************************
+"""
+
+import pandas as pd
+import matplotlib.pyplot as plt
+pd.options.mode.chained_assignment = None
+
+#%% ************ FUCNTIONS ***************************************************
+
+# windows ---------------------------------------------------------------
+def window_limit(l, window_assumptions):
+    """
+    define limit u value from which on window is retrofitted
+    """
+    m = (window_assumptions.diff()["u_limit"] /
+         window_assumptions.diff()["strength"]).dropna().iloc[0]
+    a = window_assumptions["u_limit"][0] - m * window_assumptions["strength"][0]
+    return m*l + a
+
+def u_retro_window(l, window_assumptions):
+    """
+    define retrofitting value depending on renovation strength
+    """
+    m = (window_assumptions.diff()["u_value"] /
+         window_assumptions.diff()["strength"]).dropna().iloc[0]
+    a = window_assumptions["u_value"][0] - m * window_assumptions["strength"][0]
+    return max(m*l + a, 0.8)
+
+def window_cost(u, cost_retro, window_assumptions):
+    """
+    get costs for new windows depending on u value
+
+    """
+    m = (window_assumptions.diff()["cost"] /
+         window_assumptions.diff()["u_value"]).dropna().iloc[0]
+    a = window_assumptions["cost"][0] - m * window_assumptions["u_value"][0]
+    window_cost = m*u + a
+    if annualise_cost:
+        window_cost = window_cost * interest_rate / (1 - (1 + interest_rate)
+                      ** -cost_retro.loc["Windows", "life_time"])
+    return window_cost
+
+#  functions for intermediate steps (~l, ~area)  -----------------------------
+def calculate_new_u(u_values, l, l_weight, k=0.035):
+    """
+    calculate U-values after building retrofitting, depending on the old
+    U-values (u_values).
+    They depend for the components Roof, Wall, Floor on the additional
+    insulation thickness (l), and the weighting for the corresponding
+    component (l_weight).
+    Windows are renovated to new ones with U-value (function: u_retro_window(l))
+    only if the are worse insulated than a certain limit value
+    (function: window_limit).
+
+    Parameters
+    ----------
+    u_values: pd.DataFrame
+    l: string
+    l_weight: pd.DataFrame (component, weight)
+    k: thermal conductivity
+
+    """
+    return u_values.apply(lambda x:
+                                 k / ((k / x.value) +
+                                      (float(l) * l_weight.loc[x.type][0]))
+                                 if x.type!="Windows"
+                             else (min(x.value, u_retro_window(float(l), window_assumptions))
+                                   if x.value>window_limit(float(l), window_assumptions) else x.value),
+                             axis=1)
+
+def calculate_dE(u_values, l, average_surface_w):
+    """
+    returns energy demand after retrofit (per unit of unrefurbished energy
+    demand) depending on current and retrofitted U-values, this energy demand
+    is weighted depending on the average surface of each component for the
+    building type of the assumend subsector
+    """
+    return u_values.apply(lambda x: x[l] / x.value *
+                                     average_surface_w.loc[x.assumed_subsector,
+                                                           x.type],
+                                     axis=1)
+
+
+def calculate_costs(u_values, l, cost_retro, average_surface):
+    """
+    returns costs for a given retrofitting strength weighted by the average
+    surface/volume ratio of the component for each building type
+    """
+    return u_values.apply(lambda x: (cost_retro.loc[x.type, "cost_var"] *
+                                     100 * float(l) * l_weight.loc[x.type][0]
+                                     + cost_retro.loc[x.type, "cost_fix"]) *
+                              average_surface.loc[x.assumed_subsector, x.type] /
+                              average_surface.loc[x.assumed_subsector, "surface"]
+                              if x.type!="Windows"
+                              else (window_cost(x[l], cost_retro, window_assumptions) *
+                                               average_surface.loc[x.assumed_subsector, x.type] /
+                                               average_surface.loc[x.assumed_subsector, "surface"]
+                                    if x.value>window_limit(float(l), window_assumptions) else 0),
+                             axis=1)
+
+
+# ---------------------------------------------------------------------------
+def prepare_building_stock_data():
+    """
+    reads building stock data and cleans up the format, returns
+    --------
+    u_values:          pd.DataFrame current U-values
+    average_surface:   pd.DataFrame (index= building type,
+                                     columns = [surface [m],height [m],
+                                                components area [m^2]])
+    average_surface_w: pd.DataFrame weighted share of the components per
+                       building type
+    area_tot:          heated floor area per country and sector [Mm²]
+    area:              heated floor area [Mm²] for country, sector, building
+                       type and period
+
+    """
+
+    building_data = pd.read_csv(snakemake.input.building_stock,
+                                usecols=list(range(13)))
+
+    # standardize data
+    building_data["type"].replace(
+                {'Covered area: heated  [Mm²]': 'Heated area [Mm²]',
+                 'Windows ': 'Windows',
+                 'Walls ': 'Walls',
+                 'Roof ': 'Roof',
+                 'Floor ': 'Floor'}, inplace=True)
+
+    building_data.country_code = building_data.country_code.str.upper()
+    building_data["subsector"].replace({'Hotels and Restaurants':
+                                        'Hotels and restaurants'}, inplace=True)
+    building_data["sector"].replace({'Residential sector': 'residential',
+                                     'Service sector': 'services'},
+                                    inplace=True)
+    # extract u-values
+    u_values = building_data[(building_data.feature.str.contains("U-values"))
+                             & (building_data.subsector != "Total")]
+
+    components = list(u_values.type.unique())
+
+    country_iso_dic = building_data.set_index("country")["country_code"].to_dict()
+
+    # add missing /rename countries
+    country_iso_dic.update({'Norway': 'NO',
+                            'Iceland': 'IS',
+                            'Montenegro': 'ME',
+                            'Serbia': 'RS',
+                            'Albania': 'AL',
+                            'United Kingdom': 'GB',
+                            'Bosnia and Herzegovina': 'BA',
+                            'Switzerland': 'CH'})
+
+    # heated floor area ----------------------------------------------------------
+    area = building_data[(building_data.type == 'Heated area [Mm²]') &
+                         (building_data.subsector != "Total")]
+    area_tot = area.groupby(["country", "sector"]).sum()
+    area["weight"] = area.apply(lambda x: x.value /
+                                          area_tot.value.loc[(x.country, x.sector)],
+                                axis=1)
+    area = area.groupby(['country', 'sector', 'subsector', 'bage']).sum()
+    area_tot.rename(index=country_iso_dic, inplace=True)
+
+    # add for some missing countries floor area from other data sources
+    area_missing = pd.read_csv(snakemake.input.floor_area_missing,
+                               index_col=[0, 1], usecols=[0, 1, 2, 3])
+    area_tot = area_tot.append(area_missing.unstack(level=-1).dropna().stack())
+    area_tot = area_tot.loc[~area_tot.index.duplicated(keep='last')]
+
+    # for still missing countries calculate floor area by population size
+    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()
+
+    area_per_pop = area_tot.unstack().apply(lambda x: x / ct_total[x.index])
+    missing_area_ct = ct_total.index.difference(area_tot.index.levels[0])
+    for ct in missing_area_ct:
+        averaged_data = pd.DataFrame(
+            area_per_pop.value.reindex(map_for_missings[ct]).mean()
+            * ct_total[ct],
+            columns=["value"])
+        index = pd.MultiIndex.from_product([[ct], averaged_data.index.to_list()])
+        averaged_data.index = index
+        averaged_data["estimated"] = 1
+        if ct not in area_tot.index.levels[0]:
+            area_tot = area_tot.append(averaged_data, sort=True)
+        else:
+            area_tot.loc[averaged_data.index] = averaged_data
+
+    # u_values for Poland are missing -> take them from eurostat -----------
+    u_values_PL = pd.read_csv(snakemake.input.u_values_PL)
+    area_PL = area.loc["Poland"].reset_index()
+    data_PL = pd.DataFrame(columns=u_values.columns, index=area_PL.index)
+    data_PL["country"] = "Poland"
+    data_PL["country_code"] = "PL"
+    # data from area
+    for col in ["sector", "subsector", "bage"]:
+        data_PL[col] = area_PL[col]
+    data_PL["btype"] = area_PL["subsector"]
+
+    data_PL_final = pd.DataFrame()
+    for component in components:
+        data_PL["type"] = component
+        data_PL["value"] = data_PL.apply(lambda x: u_values_PL[(u_values_PL.component==component)
+                                                               & (u_values_PL.sector==x["sector"])]
+                                         [x["bage"]].iloc[0], axis=1)
+        data_PL_final = data_PL_final.append(data_PL)
+
+    u_values = pd.concat([u_values,
+                          data_PL_final]).reset_index(drop=True)
+
+    # clean data ---------------------------------------------------------------
+    # smallest possible today u values for windows 0.8 (passive house standard)
+    # maybe the u values for the glass and not the whole window including frame
+    # for those types assumed in the dataset
+    u_values[(u_values.type=="Windows") & (u_values.value<0.8)]["value"] = 0.8
+    # drop unnecessary columns
+    u_values.drop(['topic', 'feature','detail', 'estimated','unit'],
+                  axis=1, inplace=True, errors="ignore")
+    #  only take in config.yaml specified countries into account
+    countries = ct_total.index
+    area_tot = area_tot.loc[countries]
+
+    # average component surface --------------------------------------------------
+    average_surface = (pd.read_csv(snakemake.input.average_surface,
+                                   nrows=3,
+                                   header=1,
+                                   index_col=0).rename(
+                                       {'Single/two family house': 'Single family- Terraced houses',
+                                        'Large apartment house': 'Multifamily houses',
+                                        'Apartment house': 'Appartment blocks'},
+                                    axis="index")).iloc[:, :6]
+    average_surface.columns = ["surface", "height", "Roof",
+                               "Walls", "Floor", "Windows"]
+    # get area share of component
+    average_surface_w = average_surface[components].apply(lambda x: x / x.sum(),
+                                                          axis=1)
+
+    return (u_values, average_surface,
+            average_surface_w, area_tot, area, country_iso_dic, countries)
+
+
+def prepare_cost_retro():
+    """
+    read and prepare retro costs, annualises them if annualise_cost=True
+    """
+    cost_retro = pd.read_csv(snakemake.input.cost_germany,
+                         nrows=4, index_col=0, usecols=[0, 1, 2, 3])
+    cost_retro.index = cost_retro.index.str.capitalize()
+    cost_retro.rename(index={"Window": "Windows", "Wall": "Walls"}, inplace=True)
+
+    window_assumptions = pd.read_csv(snakemake.input.window_assumptions,
+                                     skiprows=[1], usecols=[0,1,2,3], nrows=2)
+
+    if annualise_cost:
+        cost_retro[["cost_fix", "cost_var"]] = (cost_retro[["cost_fix", "cost_var"]]
+                                                .apply(lambda x: x * interest_rate /
+                                                       (1 - (1 + interest_rate)
+                                                        ** -cost_retro.loc[x.index,
+                                                                           "life_time"])))
+
+    return cost_retro, window_assumptions
+
+
+def calculate_cost_energy_curve(u_values, l_strength, l_weight, average_surface_w,
+                                average_surface, area, country_iso_dic,
+                                countries):
+    """
+    returns energy demand per unit of unrefurbished (dE) and cost for given
+    renovation strength (l_strength), data for missing countries is
+    approximated by countries with similar building stock (dict:map_for_missings)
+
+    parameter
+    -------- input -----------
+    u_values:          pd.DataFrame current U-values
+    l_strength:        list of strings (strength of retrofitting)
+    l_weight:          pd.DataFrame (component, weight)
+    average_surface:   pd.DataFrame (index= building type,
+                                     columns = [surface [m],height [m],
+                                                components area [m^2]])
+    average_surface_w: pd.DataFrame weighted share of the components per
+                       building type
+    area:              heated floor area [Mm²] for country, sector, building
+                       type and period
+    country_iso_dic:   dict (maps country name to 2-letter-iso-code)
+    countries:         pd.Index (specified countries in config.yaml)
+    -------- output ----------
+    res:               pd.DataFrame(index=pd.MultiIndex([country, sector]),
+                                    columns=pd.MultiIndex([(dE/cost), l_strength]))
+    """
+
+    energy_saved = u_values[['country', 'sector', 'subsector', 'bage', 'type']]
+    costs = u_values[['country', 'sector', 'subsector', 'bage', 'type']]
+
+    for l in l_strength:
+        u_values[l] = calculate_new_u(u_values, l, l_weight)
+        energy_saved[l] = calculate_dE(u_values, l, average_surface_w)
+        costs[l] = calculate_costs(u_values, l, cost_retro, average_surface)
+
+    # energy and costs per country, sector, subsector and year
+    e_tot = energy_saved.groupby(['country', 'sector', 'subsector', 'bage']).sum()
+    cost_tot = costs.groupby(['country', 'sector', 'subsector', 'bage']).sum()
+
+    # weighting by area -> energy and costs per country and sector
+    # in case of missing data first concat
+    energy_saved = pd.concat([e_tot, area.weight], axis=1)
+    cost_res = pd.concat([cost_tot, area.weight], axis=1)
+    energy_saved = (energy_saved.apply(lambda x: x * x.weight, axis=1)
+                    .groupby(level=[0, 1]).sum())
+    cost_res = (cost_res.apply(lambda x: x * x.weight, axis=1)
+                .groupby(level=[0, 1]).sum())
+
+    res = pd.concat([energy_saved[l_strength], cost_res[l_strength]],
+                    axis=1, keys=["dE", "cost"])
+    res.rename(index=country_iso_dic, inplace=True)
+
+    res = res.loc[countries]
+
+    # map missing countries
+    for ct in map_for_missings.keys():
+        averaged_data = pd.DataFrame(res.loc[map_for_missings[ct], :].mean(level=1))
+        index = pd.MultiIndex.from_product([[ct], averaged_data.index.to_list()])
+        averaged_data.index = index
+        if ct not in res.index.levels[0]:
+            res = res.append(averaged_data)
+        else:
+            res.loc[averaged_data.index] = averaged_data
+
+    return res
+
+
+# %% **************** MAIN ************************************************
+if __name__ == "__main__":
+    #  for testing
+    if 'snakemake' not in globals():
+        import yaml
+        import os
+        from vresutils.snakemake import MockSnakemake
+        snakemake = MockSnakemake(
+            wildcards=dict(
+                network='elec',
+                simpl='',
+                clusters='37',
+                lv='1',
+                opts='Co2L-3H',
+                sector_opts="[Co2L0p0-168H-T-H-B-I]"),
+            input=dict(
+                building_stock="data/retro/data_building_stock.csv",
+                u_values_PL="data/retro/u_values_poland.csv",
+                tax_w="data/retro/electricity_taxes_eu.csv",
+                construction_index="data/retro/comparative_level_investment.csv",
+                average_surface="data/retro/average_surface_components.csv",
+                floor_area_missing="data/retro/floor_area_missing.csv",
+                clustered_pop_layout="resources/pop_layout_{network}_s{simpl}_{clusters}.csv",
+                cost_germany="data/retro/retro_cost_germany.csv",
+                window_assumptions="data/retro/window_assumptions.csv"),
+            output=dict(
+                retro_cost="resources/retro_cost_{network}_s{simpl}_{clusters}.csv",
+                floor_area="resources/floor_area_{network}_s{simpl}_{clusters}.csv")
+        )
+        with open('config.yaml', encoding='utf8') as f:
+            snakemake.config = yaml.safe_load(f)
+
+#  ******** (1) ASSUMPTIONS - PARAMETERS **********************************
+    retro_opts =  snakemake.config["sector"]["retrofitting"]
+    interest_rate = retro_opts["interest_rate"]
+    annualise_cost = retro_opts["annualise_cost"]  # annualise the investment costs
+    tax_weighting = retro_opts["tax_weighting"]   # weight costs depending on taxes in countries
+    construction_index = retro_opts["construction_index"]   # weight costs depending on labour/material costs per ct
+    # additional insulation thickness, determines maximum possible savings
+    l_strength = retro_opts["l_strength"]
+
+    k = 0.035   # thermal conductivity standard value
+    # strenght of relative retrofitting depending on the component
+    # determined by historical data of insulation thickness for retrofitting
+    l_weight = pd.DataFrame({"weight": [1.95, 1.48, 1.]},
+                            index=["Roof", "Walls", "Floor"])
+
+
+    # mapping missing countries by neighbours
+    map_for_missings = {
+        "AL": ["BG", "RO", "GR"],
+        "BA": ["HR"],
+        "RS": ["BG", "RO", "HR", "HU"],
+        "MK": ["BG", "GR"],
+        "ME": ["BA", "AL", "RS", "HR"],
+        "CH": ["SE", "DE"],
+        "NO": ["SE"],
+        }
+
+# %% ************ (2) DATA ***************************************************
+
+    # building stock data -----------------------------------------------------
+    (u_values, average_surface, average_surface_w,
+     area_tot, area, country_iso_dic, countries) = prepare_building_stock_data()
+
+    # costs for retrofitting -------------------------------------------------
+    cost_retro, window_assumptions = prepare_cost_retro()
+
+    # weightings of costs
+    if construction_index:
+        cost_w = pd.read_csv(snakemake.input.construction_index,
+                             skiprows=3, nrows=32, index_col=0)
+        # since German retrofitting costs are assumed
+        cost_w = ((cost_w["2018"] / cost_w.loc["Germany", "2018"])
+                  .rename(index=country_iso_dic))
+
+    if tax_weighting:
+        tax_w = pd.read_csv(snakemake.input.tax_w,
+                            header=12, nrows=39, index_col=0, usecols=[0, 4])
+        tax_w.rename(index=country_iso_dic, inplace=True)
+        tax_w = tax_w.apply(pd.to_numeric, errors='coerce').iloc[:, 0]
+        tax_w.dropna(inplace=True)
+
+# %% ********** (3) CALCULATE COST-ENERGY-CURVES ****************************
+
+    # for missing weighting of surfaces of building types assume MultiFamily houses
+    u_values["assumed_subsector"] = u_values.subsector
+    u_values.assumed_subsector[
+        ~u_values.subsector.isin(average_surface.index)] = 'Multifamily houses'
+
+    dE_and_cost =  calculate_cost_energy_curve(u_values, l_strength, l_weight,
+                                          average_surface_w, average_surface, area,
+                                          country_iso_dic, countries)
+
+    # weights costs after construction index
+    if construction_index:
+        for ct in list(map_for_missings.keys() - cost_w.index):
+            cost_w.loc[ct] = cost_w.reindex(index=map_for_missings[ct]).mean()
+        dE_and_cost.cost = dE_and_cost.cost.apply(lambda x: x * cost_w[x.index.levels[0]])
+
+    # weights cost depending on country taxes
+    if tax_weighting:
+        for ct in list(map_for_missings.keys() - tax_w.index):
+            tax_w[ct] = tax_w.reindex(index=map_for_missings[ct]).mean()
+        dE_and_cost.cost = dE_and_cost.cost.apply(lambda x: x * tax_w[x.index.levels[0]])
+
+    # get share of residential and sevice floor area
+    sec_w = (area_tot / area_tot.groupby(["country"]).sum())["value"]
+    # get the total cost-energy-savings weight by sector area
+    tot = dE_and_cost.apply(lambda col: col * sec_w, axis=0).groupby(level=0).sum()
+    tot.set_index(pd.MultiIndex.from_product([list(tot.index), ["tot"]]),
+                  inplace=True)
+    dE_and_cost = dE_and_cost.append(tot).unstack().stack()
+
+    summed_area = pd.DataFrame(area_tot.groupby("country").sum())
+    summed_area.set_index(pd.MultiIndex.from_product(
+                          [list(summed_area.index), ["tot"]]), inplace=True)
+    area_tot = area_tot.append(summed_area).unstack().stack()
+
+# %% ******* (4) SAVE   ************************************************
+
+    dE_and_cost.to_csv(snakemake.output.retro_cost)
+    area_tot.to_csv(snakemake.output.floor_area)
+
+
+

scripts/make_summary.py

@@ -592,7 +592,8 @@ if __name__ == "__main__":
     costs_db = prepare_costs(snakemake.input.costs,
                              snakemake.config['costs']['USD2013_to_EUR2013'],
                              snakemake.config['costs']['discountrate'],
-                             Nyears)
+                             Nyears,
+                             snakemake.config['costs']['lifetime'])
 
     df = make_summaries(networks_dict)
 

scripts/prepare_sector_network.py

@@ -114,6 +114,13 @@ def update_wind_solar_costs(n,costs):
             n.generators.loc[n.generators.carrier==tech,'capital_cost'] = capital_cost.rename(index=lambda node: node + ' ' + tech)
 
 
+def retro_exogen(demand, dE):
+    """
+    reduces space heat demand exogenously
+    demand: current space heat demand
+    dE: energy savings
+    """
+    return demand * (1-dE)
 def add_carrier_buses(n, carriers):
     """
     Add buses to connect e.g. coal, nuclear and oil plants
@@ -451,8 +458,8 @@ def prepare_data(network):
 
     ## Get overall demand curve for all vehicles
 
-    dir_name = "data/emobility/"
+    traffic = pd.read_csv(snakemake.input.traffic_data + "KFZ__count",
-    traffic = pd.read_csv(os.path.join(dir_name,"KFZ__count"),skiprows=2)["count"]
+                          skiprows=2)["count"]
 
     #Generate profiles
     transport_shape = generate_periodic_profiles(dt_index=network.snapshots.tz_localize("UTC"),
@@ -506,7 +513,8 @@ def prepare_data(network):
 
     ## derive plugged-in availability for PKW's (cars)
 
-    traffic = pd.read_csv(os.path.join(dir_name,"Pkw__count"),skiprows=2)["count"]
+    traffic = pd.read_csv(snakemake.input.traffic_data + "Pkw__count",
+                          skiprows=2)["count"]
 
     avail_max = 0.95
 
@@ -540,7 +548,7 @@ def prepare_data(network):
 
 
 
-def prepare_costs(cost_file, USD_to_EUR, discount_rate, Nyears):
+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()
@@ -558,7 +566,7 @@ def prepare_costs(cost_file, USD_to_EUR, discount_rate, Nyears):
                           "efficiency" : 1,
                           "fuel" : 0,
                           "investment" : 0,
-                          "lifetime" : 25
+                          "lifetime" : lifetime
     })
 
     costs["fixed"] = [(annuity(v["lifetime"],v["discount rate"])+v["FOM"]/100.)*v["investment"]*Nyears for i,v in costs.iterrows()]
@@ -1077,7 +1085,21 @@ def add_heat(network):
 
     urban_fraction = options['central_fraction']*pop_layout["urban"]/(pop_layout[["urban","rural"]].sum(axis=1))
 
-    for name in ["residential rural","services rural","residential urban decentral","services urban decentral","urban central"]:
+    # building retrofitting, exogenously reduce space heat demand
+    if options["retrofitting"]["retro_exogen"]:
+        dE = options["retrofitting"]["dE"]
+        if snakemake.config["foresight"]=='myopic':
+            year = int(snakemake.wildcards.planning_horizons[-4:])
+            dE = dE[snakemake.config["scenario"]["planning_horizons"].index(year)]
+        print("retrofitting exogenously, assumed space heat reduction of ",
+              dE)
+        for sector in sectors:
+            heat_demand[sector + " space"] = heat_demand[sector + " space"].apply(lambda x: retro_exogen(x, dE))
+
+    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"
 
@@ -1100,6 +1122,7 @@ def add_heat(network):
             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']))
 
@@ -1291,64 +1314,100 @@ def add_heat(network):
                              lifetime=costs.at['micro CHP','lifetime'])
 
 
-    #NB: this currently doesn't work for pypsa-eur model
+    if options['retrofitting']['retro_endogen']:
-    if options['retrofitting']:
 
-        retro_nodes = pd.Index(["DE"])
+        print("adding retrofitting endogenously")
 
-        space_heat_demand = space_heat_demand[retro_nodes]
+        # 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()
 
-        square_metres = population[retro_nodes]/population['DE']*5.7e9   #HPI 3.4e9m^2 for DE res, 2.3e9m^2 for tert https://doi.org/10.1016/j.rser.2013.09.012
+        # retrofitting data 'retro_data' with 'costs' [EUR/m^2] and heat
-
+        # demand 'dE' [per unit of original heat demand] for each country and
-        space_peak = space_heat_demand.max()
+        # different retrofitting strengths [additional insulation thickness in m]
-
+        retro_data = pd.read_csv(snakemake.input.retro_cost_energy,
-        space_pu = space_heat_demand.divide(space_peak)
+                                 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")
 
-        network.madd('Generator',
+        # share of space heat demand 'w_space' of total heat demand
-                     retro_nodes,
+        w_space = {}
-                     suffix=' retrofitting I',
+        for sector in sectors:
-                     bus=retro_nodes+' heat',
+            w_space[sector] = heat_demand_r[sector + " space"] / \
-                     carrier="retrofitting",
+                (heat_demand_r[sector + " space"] + heat_demand_r[sector + " water"])
-                     p_nom_extendable=True,
+        w_space["tot"] = ((heat_demand_r["services space"] +
-                     p_nom_max=options['retroI-fraction']*space_peak*(1-urban_fraction),
+                           heat_demand_r["residential space"]) /
-                     p_max_pu=space_pu,
+                           heat_demand_r.groupby(level=[1], axis=1).sum())
-                     p_min_pu=space_pu,
-                     capital_cost=options['retrofitting-cost_factor']*costs.at['retrofitting I','fixed']*square_metres/(options['retroI-fraction']*space_peak))
 
-        network.madd('Generator',
-                     retro_nodes,
-                     suffix=' retrofitting II',
-                     bus=retro_nodes+' heat',
-                     carrier="retrofitting",
-                     p_nom_extendable=True,
-                     p_nom_max=options['retroII-fraction']*space_peak*(1-urban_fraction),
-                     p_max_pu=space_pu,
-                     p_min_pu=space_pu,
-                     capital_cost=options['retrofitting-cost_factor']*costs.at['retrofitting II','fixed']*square_metres/(options['retroII-fraction']*space_peak))
 
-        network.madd('Generator',
+        for name in network.loads[network.loads.carrier.isin([x + " heat" for x in heat_systems])].index:
-                     retro_nodes,
+
-                     suffix=' urban retrofitting I',
+            node = network.buses.loc[name, "location"]
-                     bus=retro_nodes+' urban heat',
+            ct = pop_layout.loc[node, "ct"]
-                     carrier="retrofitting",
+
-                     p_nom_extendable=True,
+            # weighting 'f' depending on the size of the population at the node
-                     p_nom_max=options['retroI-fraction']*space_peak*urban_fraction,
+            f = urban_fraction[node] if "urban" in name else (1-urban_fraction[node])
-                     p_max_pu=space_pu,
+            if f == 0:
-                     p_min_pu=space_pu,
+                continue
-                     capital_cost=options['retrofitting-cost_factor']*costs.at['retrofitting I','fixed']*square_metres/(options['retroI-fraction']*space_peak))
+            # 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'])
 
-        network.madd('Generator',
-                     retro_nodes,
-                     suffix=' urban retrofitting II',
-                     bus=retro_nodes+' urban heat',
-                     carrier="retrofitting",
-                     p_nom_extendable=True,
-                     p_nom_max=options['retroII-fraction']*space_peak*urban_fraction,
-                     p_max_pu=space_pu,
-                     p_min_pu=space_pu,
-                     capital_cost=options['retrofitting-cost_factor']*costs.at['retrofitting II','fixed']*square_metres/(options['retroII-fraction']*space_peak))
 
 
 def create_nodes_for_heat_sector():
@@ -1782,25 +1841,44 @@ if __name__ == "__main__":
             wildcards=dict(network='elec', simpl='', clusters='37', lv='1.0',
                            opts='', planning_horizons='2020',
                            sector_opts='Co2L0-168H-T-H-B-I-solar3-dist1'),
-            input=dict(network='pypsa-eur/networks/{network}_s{simpl}_{clusters}_ec_lv{lv}_{opts}.nc',
+            input=dict(network='../pypsa-eur/networks/{network}_s{simpl}_{clusters}_ec_lv{lv}_{opts}.nc',
-                       timezone_mappings='pypsa-eur-sec/data/timezone_mappings.csv',
+                       energy_totals_name='resources/energy_totals.csv',
-                       clustered_pop_layout='pypsa-eur-sec/resources/pop_layout_{network}_s{simpl}_{clusters}.csv',
+                       co2_totals_name='resources/co2_totals.csv',
-                       costs='technology-data/outputs/costs_{planning_horizons}.csv',
+                       transport_name='resources/transport_data.csv',
-                       profile_offwind_ac='pypsa-eur/resources/profile_offwind-ac.nc',
+                       biomass_potentials='resources/biomass_potentials.csv',
-                       profile_offwind_dc='pypsa-eur/resources/profile_offwind-dc.nc',
+                       biomass_transport='data/biomass/biomass_transport_costs.csv',
-                       busmap_s='pypsa-eur/resources/busmap_{network}_s{simpl}.csv',
+                       timezone_mappings='data/timezone_mappings.csv',
-                       busmap='pypsa-eur/resources/busmap_{network}_s{simpl}_{clusters}.csv',
+                       heat_profile="data/heat_load_profile_BDEW.csv",
-                       cop_air_total='pypsa-eur-sec/resources/cop_air_total_{network}_s{simpl}_{clusters}.nc',
+                       costs="../technology-data/outputs/costs_{planning_horizons}.csv",
-                       cop_soil_total='pypsa-eur-sec/resources/cop_soil_total_{network}_s{simpl}_{clusters}.nc',
+	               h2_cavern = "data/hydrogen_salt_cavern_potentials.csv",
-                       solar_thermal_total='pypsa-eur-sec/resources/solar_thermal_total_{network}_s{simpl}_{clusters}.nc',
+                       profile_offwind_ac="../pypsa-eur/resources/profile_offwind-ac.nc",
-                       energy_totals_name='pypsa-eur-sec/data/energy_totals.csv',
+                       profile_offwind_dc="../pypsa-eur/resources/profile_offwind-dc.nc",
-                       heat_demand_total='pypsa-eur-sec/resources/heat_demand_total_{network}_s{simpl}_{clusters}.nc',
+                       clustermaps='../pypsa-eur/resources/clustermaps_{network}_s{simpl}_{clusters}.h5',
-                       heat_profile='pypsa-eur-sec/data/heat_load_profile_BDEW.csv',
+                       clustered_pop_layout="resources/pop_layout_{network}_s{simpl}_{clusters}.csv",
-                       transport_name='pypsa-eur-sec/data/transport_data.csv',
+                       simplified_pop_layout="resources/pop_layout_{network}_s{simpl}.csv",
-                       temp_air_total='pypsa-eur-sec/resources/temp_air_total_{network}_s{simpl}_{clusters}.nc',
+                       industrial_demand="resources/industrial_energy_demand_{network}_s{simpl}_{clusters}.csv",
-                       co2_totals_name='pypsa-eur-sec/data/co2_totals.csv',
+                       heat_demand_urban="resources/heat_demand_urban_{network}_s{simpl}_{clusters}.nc",
-                       biomass_potentials='pypsa-eur-sec/data/biomass_potentials.csv',
+                       heat_demand_rural="resources/heat_demand_rural_{network}_s{simpl}_{clusters}.nc",
-                       industrial_demand='pypsa-eur-sec/resources/industrial_demand_{network}_s{simpl}_{clusters}.csv',),
+                       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",
+                       traffic_data = "data/emobility/",
+                       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=['pypsa-eur-sec/results/test/prenetworks/{network}_s{simpl}_{clusters}_lv{lv}__{sector_opts}_{co2_budget_name}_{planning_horizons}.nc']
         )
         import yaml
@@ -1834,7 +1912,8 @@ if __name__ == "__main__":
     costs = prepare_costs(snakemake.input.costs,
                           snakemake.config['costs']['USD2013_to_EUR2013'],
                           snakemake.config['costs']['discountrate'],
-                          Nyears)
+                          Nyears,
+                          snakemake.config['costs']['lifetime'])
 
     remove_elec_base_techs(n)