From cfe44ce702331a7dd749cf28e41f1d9f10eb9b7f Mon Sep 17 00:00:00 2001 From: Tom Brown Date: Thu, 18 Apr 2019 11:39:17 +0200 Subject: [PATCH] Copy solve_network.py and *_summary.py from pypsa-eur sector branch --- scripts/make_summary.py | 445 +++++++++++++++++++++++++++++++++++++++ scripts/plot_summary.py | 187 ++++++++++++++++ scripts/solve_network.py | 375 +++++++++++++++++++++++++++++++++ 3 files changed, 1007 insertions(+) create mode 100644 scripts/make_summary.py create mode 100644 scripts/plot_summary.py create mode 100644 scripts/solve_network.py diff --git a/scripts/make_summary.py b/scripts/make_summary.py new file mode 100644 index 00000000..4a7b348e --- /dev/null +++ b/scripts/make_summary.py @@ -0,0 +1,445 @@ + +from six import iteritems + +import sys + +sys.path = ['/home/vres/lib/python3.5/site-packages'] + sys.path + +import pandas as pd + +import numpy as np + +import pypsa + +from vresutils.costdata import annuity + +from prepare_network import generate_periodic_profiles + +from add_electricity import load_costs + +import yaml + +idx = pd.IndexSlice + +opt_name = {"Store": "e", "Line" : "s", "Transformer" : "s"} + +#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["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["p2"] = ["series","MW",0.,"2nd bus output","Output"] +override_component_attrs["Link"].loc["p3"] = ["series","MW",0.,"3rd bus output","Output"] + + + + +def assign_carriers(n): + + if "carrier" not in n.loads: + n.loads["carrier"] = "electricity" + for carrier in ["transport","heat","urban heat"]: + n.loads.loc[n.loads.index.str.contains(carrier),"carrier"] = carrier + + if "carrier" not in n.lines: + n.lines["carrier"] = "AC" + + if n.stores.loc["EU gas Store","carrier"] == "": + n.stores.loc["EU gas Store","carrier"] = "gas Store" + + +def calculate_costs(n,label,costs): + + for c in n.iterate_components(n.branch_components|n.controllable_one_port_components^{"Load"}): + capital_costs = c.df.capital_cost*c.df[opt_name.get(c.name,"p") + "_nom_opt"] + capital_costs_grouped = capital_costs.groupby(c.df.carrier).sum() + + costs = costs.reindex(costs.index|pd.MultiIndex.from_product([[c.list_name],["capital"],capital_costs_grouped.index])) + + costs.loc[idx[c.list_name,"capital",list(capital_costs_grouped.index)],label] = capital_costs_grouped.values + + if c.name == "Link": + p = c.pnl.p0.multiply(n.snapshot_weightings,axis=0).sum() + elif c.name == "Line": + continue + elif c.name == "StorageUnit": + p_all = c.pnl.p.multiply(n.snapshot_weightings,axis=0) + p_all[p_all < 0.] = 0. + p = p_all.sum() + else: + p = c.pnl.p.multiply(n.snapshot_weightings,axis=0).sum() + + marginal_costs = p*c.df.marginal_cost + + marginal_costs_grouped = marginal_costs.groupby(c.df.carrier).sum() + + costs = costs.reindex(costs.index|pd.MultiIndex.from_product([[c.list_name],["marginal"],marginal_costs_grouped.index])) + + costs.loc[idx[c.list_name,"marginal",list(marginal_costs_grouped.index)],label] = marginal_costs_grouped.values + + #add back in costs of links if there is a line volume limit + if label[1] != "opt": + costs.loc[("links-added","capital","transmission lines"),label] = ((costs_db.at['HVDC overhead', 'capital_cost']*n.links.length + costs_db.at['HVDC inverter pair', 'capital_cost'])*n.links.p_nom_opt)[n.links.carrier == "DC"].sum() + costs.loc[("lines-added","capital","transmission lines"),label] = costs_db.at["HVAC overhead", "capital_cost"]*(n.lines.length*n.lines.s_nom_opt).sum() + else: + costs.loc[("links-added","capital","transmission lines"),label] = (costs_db.at['HVDC inverter pair', 'capital_cost']*n.links.p_nom_opt)[n.links.carrier == "DC"].sum() + + + #add back in all hydro + #costs.loc[("storage_units","capital","hydro"),label] = (0.01)*2e6*n.storage_units.loc[n.storage_units.group=="hydro","p_nom"].sum() + #costs.loc[("storage_units","capital","PHS"),label] = (0.01)*2e6*n.storage_units.loc[n.storage_units.group=="PHS","p_nom"].sum() + #costs.loc[("generators","capital","ror"),label] = (0.02)*3e6*n.generators.loc[n.generators.group=="ror","p_nom"].sum() + + return costs + + + +def calculate_curtailment(n,label,curtailment): + + avail = n.generators_t.p_max_pu.multiply(n.generators.p_nom_opt).sum().groupby(n.generators.carrier).sum() + used = n.generators_t.p.sum().groupby(n.generators.carrier).sum() + + curtailment[label] = (((avail - used)/avail)*100).round(3) + + return curtailment + +def calculate_energy(n,label,energy): + + for c in n.iterate_components(n.one_port_components|n.branch_components): + + if c.name in n.one_port_components: + c_energies = c.pnl.p.multiply(n.snapshot_weightings,axis=0).sum().multiply(c.df.sign).groupby(c.df.carrier).sum() + else: + c_energies = (-c.pnl.p1.multiply(n.snapshot_weightings,axis=0).sum() - c.pnl.p0.multiply(n.snapshot_weightings,axis=0).sum()).groupby(c.df.carrier).sum() + + energy = energy.reindex(energy.index|pd.MultiIndex.from_product([[c.list_name],c_energies.index])) + + energy.loc[idx[c.list_name,list(c_energies.index)],label] = c_energies.values + + return energy + + +def calculate_supply(n,label,supply): + """calculate the max dispatch of each component at the buses where the loads are attached""" + + load_types = n.loads.carrier.value_counts().index + + for i in load_types: + + buses = n.loads.bus[n.loads.carrier == i].values + + bus_map = pd.Series(False,index=n.buses.index) + + bus_map.loc[buses] = True + + for c in n.iterate_components(n.one_port_components): + + items = c.df.index[c.df.bus.map(bus_map)] + + if len(items) == 0: + continue + + s = c.pnl.p[items].max().multiply(c.df.loc[items,'sign']).groupby(c.df.loc[items,'carrier']).sum() + + supply = supply.reindex(supply.index|pd.MultiIndex.from_product([[i],[c.list_name],s.index])) + supply.loc[idx[i,c.list_name,list(s.index)],label] = s.values + + + for c in n.iterate_components(n.branch_components): + + for end in ["0","1"]: + + items = c.df.index[c.df["bus" + end].map(bus_map)] + + if len(items) == 0: + continue + + #lots of sign compensation for direction and to do maximums + s = (-1)**(1-int(end))*((-1)**int(end)*c.pnl["p"+end][items]).max().groupby(c.df.loc[items,'carrier']).sum() + + supply = supply.reindex(supply.index|pd.MultiIndex.from_product([[i],[c.list_name],s.index])) + supply.loc[idx[i,c.list_name,list(s.index)],label] = s.values + + return supply + +def calculate_supply_energy(n,label,supply_energy): + """calculate the total dispatch of each component at the buses where the loads are attached""" + + load_types = n.loads.carrier.value_counts().index + + for i in load_types: + + buses = n.loads.bus[n.loads.carrier == i].values + + bus_map = pd.Series(False,index=n.buses.index) + + bus_map.loc[buses] = True + + for c in n.iterate_components(n.one_port_components): + + items = c.df.index[c.df.bus.map(bus_map)] + + if len(items) == 0: + continue + + s = c.pnl.p[items].sum().multiply(c.df.loc[items,'sign']).groupby(c.df.loc[items,'carrier']).sum() + + supply_energy = supply_energy.reindex(supply_energy.index|pd.MultiIndex.from_product([[i],[c.list_name],s.index])) + supply_energy.loc[idx[i,c.list_name,list(s.index)],label] = s.values + + + for c in n.iterate_components(n.branch_components): + + for end in ["0","1"]: + + items = c.df.index[c.df["bus" + end].map(bus_map)] + + if len(items) == 0: + continue + + s = (-1)*c.pnl["p"+end][items].sum().groupby(c.df.loc[items,'carrier']).sum() + + supply_energy = supply_energy.reindex(supply_energy.index|pd.MultiIndex.from_product([[i],[c.list_name],s.index])) + supply_energy.loc[idx[i,c.list_name,list(s.index)],label] = s.values + + return supply_energy + +def calculate_metrics(n,label,metrics): + + metrics = metrics.reindex(metrics.index|pd.Index(["line_volume","line_volume_limit","line_volume_AC","line_volume_DC","line_volume_shadow","co2_shadow"])) + + metrics.at["line_volume_DC",label] = (n.links.length*n.links.p_nom_opt)[n.links.carrier == "DC"].sum() + metrics.at["line_volume_AC",label] = (n.lines.length*n.lines.s_nom_opt).sum() + metrics.at["line_volume",label] = metrics.loc[["line_volume_AC","line_volume_DC"],label].sum() + + if hasattr(n,"line_volume_limit"): + metrics.at["line_volume_limit",label] = n.line_volume_limit + metrics.at["line_volume_shadow",label] = n.line_volume_limit_dual + + if "CO2Limit" in n.global_constraints.index: + metrics.at["co2_shadow",label] = n.global_constraints.at["CO2Limit","mu"] + + return metrics + + +def calculate_prices(n,label,prices): + + bus_type = pd.Series(n.buses.index.str[3:],n.buses.index).replace("","electricity") + + prices = prices.reindex(prices.index|bus_type.value_counts().index) + + #WARNING: this is time-averaged, should really be load-weighted average + prices[label] = n.buses_t.marginal_price.mean().groupby(bus_type).mean() + + return prices + + + +def calculate_weighted_prices(n,label,weighted_prices): + # Warning: doesn't include storage units as loads + + + weighted_prices = weighted_prices.reindex(pd.Index(["electricity","heat","space heat","urban heat","space urban heat","gas","H2"])) + + link_loads = {"electricity" : ["heat pump", "resistive heater", "battery charger", "H2 Electrolysis"], + "heat" : ["water tanks charger"], + "urban heat" : ["water tanks charger"], + "space heat" : [], + "space urban heat" : [], + "gas" : ["OCGT","gas boiler","CHP electric","CHP heat"], + "H2" : ["Sabatier", "H2 Fuel Cell"]} + + for carrier in link_loads: + + if carrier == "electricity": + suffix = "" + elif carrier[:5] == "space": + suffix = carrier[5:] + else: + suffix = " " + carrier + + buses = n.buses.index[n.buses.index.str[2:] == suffix] + + if buses.empty: + continue + + if carrier in ["H2","gas"]: + load = pd.DataFrame(index=n.snapshots,columns=buses,data=0.) + elif carrier[:5] == "space": + load = heat_demand_df[buses.str[:2]].rename(columns=lambda i: str(i)+suffix) + else: + load = n.loads_t.p_set[buses] + + + for tech in link_loads[carrier]: + + names = n.links.index[n.links.index.to_series().str[-len(tech):] == tech] + + if names.empty: + continue + + load += n.links_t.p0[names].groupby(n.links.loc[names,"bus0"],axis=1).sum(axis=1) + + #Add H2 Store when charging + if carrier == "H2": + stores = n.stores_t.p[buses+ " Store"].groupby(n.stores.loc[buses+ " Store","bus"],axis=1).sum(axis=1) + stores[stores > 0.] = 0. + load += -stores + + weighted_prices.loc[carrier,label] = (load*n.buses_t.marginal_price[buses]).sum().sum()/load.sum().sum() + + if carrier[:5] == "space": + print(load*n.buses_t.marginal_price[buses]) + + return weighted_prices + + + + +def calculate_market_values(n, label, market_values): + # Warning: doesn't include storage units + + carrier = "AC" + + buses = n.buses.index[n.buses.carrier == carrier] + + ## First do market value of generators ## + + generators = n.generators.index[n.buses.loc[n.generators.bus,"carrier"] == carrier] + + techs = n.generators.loc[generators,"carrier"].value_counts().index + + market_values = market_values.reindex(market_values.index | techs) + + + for tech in techs: + gens = generators[n.generators.loc[generators,"carrier"] == tech] + + dispatch = n.generators_t.p[gens].groupby(n.generators.loc[gens,"bus"],axis=1).sum().reindex(columns=buses,fill_value=0.) + + revenue = dispatch*n.buses_t.marginal_price[buses] + + market_values.at[tech,label] = revenue.sum().sum()/dispatch.sum().sum() + + + + ## Now do market value of links ## + + for i in ["0","1"]: + all_links = n.links.index[n.buses.loc[n.links["bus"+i],"carrier"] == carrier] + + techs = n.links.loc[all_links,"carrier"].value_counts().index + + market_values = market_values.reindex(market_values.index | techs) + + for tech in techs: + links = all_links[n.links.loc[all_links,"carrier"] == tech] + + dispatch = n.links_t["p"+i][links].groupby(n.links.loc[links,"bus"+i],axis=1).sum().reindex(columns=buses,fill_value=0.) + + revenue = dispatch*n.buses_t.marginal_price[buses] + + market_values.at[tech,label] = revenue.sum().sum()/dispatch.sum().sum() + + return market_values + + +def calculate_price_statistics(n, label, price_statistics): + + + price_statistics = price_statistics.reindex(price_statistics.index|pd.Index(["zero_hours","mean","standard_deviation"])) + + buses = n.buses.index[n.buses.carrier == "AC"] + + threshold = 0.1 #higher than phoney marginal_cost of wind/solar + + df = pd.DataFrame(data=0.,columns=buses,index=n.snapshots) + + df[n.buses_t.marginal_price[buses] < threshold] = 1. + + price_statistics.at["zero_hours", label] = df.sum().sum()/(df.shape[0]*df.shape[1]) + + price_statistics.at["mean", label] = n.buses_t.marginal_price[buses].unstack().mean() + + price_statistics.at["standard_deviation", label] = n.buses_t.marginal_price[buses].unstack().std() + + return price_statistics + + +outputs = ["costs", + "curtailment", + "energy", + "supply", + "supply_energy", + "prices", + "weighted_prices", + "price_statistics", + "market_values", + "metrics", + ] + +def make_summaries(networks_dict): + + columns = pd.MultiIndex.from_tuples(networks_dict.keys(),names=["cluster","lv","opt"]) + + df = {} + + for output in outputs: + df[output] = pd.DataFrame(columns=columns,dtype=float) + + for label, filename in iteritems(networks_dict): + print(label, filename) + + n = pypsa.Network(filename, + override_component_attrs=override_component_attrs) + + + assign_carriers(n) + + for output in outputs: + df[output] = globals()["calculate_" + output](n, label, df[output]) + + return df + + +def to_csv(df): + + for key in df: + df[key].to_csv(snakemake.output[key]) + + +if __name__ == "__main__": + # Detect running outside of snakemake and mock snakemake for testing + if 'snakemake' not in globals(): + from vresutils import Dict + import yaml + snakemake = Dict() + with open('config.yaml') as f: + snakemake.config = yaml.load(f) + snakemake.input = Dict() + snakemake.input['heat_demand_name'] = 'data/heating/daily_heat_demand.h5' + snakemake.output = Dict() + for item in outputs: + snakemake.output[item] = snakemake.config['summary_dir'] + '/{name}/csvs/{item}.csv'.format(name=snakemake.config['run'],item=item) + + networks_dict = {(cluster,lv,opt) : + snakemake.config['results_dir'] + snakemake.config['run'] + '/postnetworks/elec_s_{cluster}_lv{lv}_{opt}.nc'\ + .format(cluster=cluster, + opt=opt, + lv=lv)\ + for cluster in snakemake.config['scenario']['clusters'] \ + for opt in snakemake.config['scenario']['opts'] \ + for lv in snakemake.config['scenario']['lv']} + print(networks_dict) + + costs_db = load_costs(Nyears=1.,tech_costs="data/costs.csv",config=snakemake.config["costs"],elec_config=snakemake.config['electricity']) + + df = make_summaries(networks_dict) + + to_csv(df) diff --git a/scripts/plot_summary.py b/scripts/plot_summary.py new file mode 100644 index 00000000..52b76141 --- /dev/null +++ b/scripts/plot_summary.py @@ -0,0 +1,187 @@ + + + +import pandas as pd + +#allow plotting without Xwindows +import matplotlib +matplotlib.use('Agg') + +import matplotlib.pyplot as plt + + + +#consolidate and rename +def rename_techs(label): + if label[:8] == "central ": + label = label[8:] + if label[:6] == "urban ": + label = label[6:] + + if "retrofitting" in label: + label = "building retrofitting" + if "H2" in label: + label = "hydrogen storage" + if "CHP" in label: + label = "CHP" + if "water tank" in label: + label = "water tanks" + if label=="water tanks": + label = "hot water storage" + if "gas" in label and label not in ["gas boiler","biogas"]: + label = "natural gas" + if "solar thermal" in label: + label = "solar thermal" + if label == "solar": + label = "solar PV" + if label == "heat pump": + label = "air heat pump" + if label == "Sabatier": + label = "methanation" + if label == "offwind": + label = "offshore wind" + if label == "onwind": + label = "onshore wind" + if label == "ror": + label = "hydroelectricity" + if label == "hydro": + label = "hydroelectricity" + if label == "PHS": + label = "hydroelectricity" + if label == "co2 Store": + label = "DAC" + if label == "co2 stored": + label = "CO2 sequestration" + if "battery" in label: + label = "battery storage" + if label in ["AC","DC","B2B"]: + label = "transmission lines" + return label + + +preferred_order = pd.Index(["transmission lines","hydroelectricity","hydro reservoir","run of river","pumped hydro storage","solid biomass","biogas","onshore wind","offshore wind","solar PV","solar thermal","solar","building retrofitting","ground heat pump","air heat pump","heat pump","resistive heater","power-to-heat","gas-to-power/heat","CHP","OCGT","gas boiler","gas","natural gas","helmeth","methanation","hydrogen storage","power-to-gas","power-to-liquid","battery storage","hot water storage","CO2 sequestration"]) + +def plot_costs(): + + + cost_df = pd.read_csv(snakemake.input.costs,index_col=list(range(3)),header=[0,1,2]) + + + df = cost_df.groupby(cost_df.index.get_level_values(2)).sum() + + #convert to billions + df = df/1e9 + + df = df.groupby(df.index.map(rename_techs)).sum() + + to_drop = df.index[df.max(axis=1) < snakemake.config['plotting']['costs_threshold']] + + print("dropping") + + print(df.loc[to_drop]) + + df = df.drop(to_drop) + + print(df.sum()) + + new_index = (preferred_order&df.index).append(df.index.difference(preferred_order)) + + new_columns = df.sum().sort_values().index + + fig, ax = plt.subplots() + fig.set_size_inches((12,8)) + + df.loc[new_index,new_columns].T.plot(kind="bar",ax=ax,stacked=True,color=[snakemake.config['plotting']['tech_colors'][i] for i in new_index]) + + + handles,labels = ax.get_legend_handles_labels() + + handles.reverse() + labels.reverse() + + ax.set_ylim([0,snakemake.config['plotting']['costs_max']]) + + ax.set_ylabel("System Cost [EUR billion per year]") + + ax.set_xlabel("") + + ax.grid(axis="y") + + ax.legend(handles,labels,ncol=4,loc="upper left") + + + fig.tight_layout() + + fig.savefig(snakemake.output.costs,transparent=True) + + +def plot_energy(): + + energy_df = pd.read_csv(snakemake.input.energy,index_col=list(range(2)),header=[0,1,2]) + + df = energy_df.groupby(energy_df.index.get_level_values(1)).sum() + + #convert MWh to TWh + df = df/1e6 + + df = df.groupby(df.index.map(rename_techs)).sum() + + to_drop = df.index[df.abs().max(axis=1) < snakemake.config['plotting']['energy_threshold']] + + print("dropping") + + print(df.loc[to_drop]) + + df = df.drop(to_drop) + + print(df.sum()) + + new_index = (preferred_order&df.index).append(df.index.difference(preferred_order)) + + new_columns = df.columns.sort_values() + + fig, ax = plt.subplots() + fig.set_size_inches((12,8)) + + df.loc[new_index,new_columns].T.plot(kind="bar",ax=ax,stacked=True,color=[snakemake.config['plotting']['tech_colors'][i] for i in new_index]) + + + handles,labels = ax.get_legend_handles_labels() + + handles.reverse() + labels.reverse() + + ax.set_ylim([snakemake.config['plotting']['energy_min'],snakemake.config['plotting']['energy_max']]) + + ax.set_ylabel("Energy [TWh/a]") + + ax.set_xlabel("") + + ax.grid(axis="y") + + ax.legend(handles,labels,ncol=4,loc="upper left") + + + fig.tight_layout() + + fig.savefig(snakemake.output.energy,transparent=True) + + +if __name__ == "__main__": + # Detect running outside of snakemake and mock snakemake for testing + if 'snakemake' not in globals(): + from vresutils import Dict + import yaml + snakemake = Dict() + with open('config.yaml') as f: + snakemake.config = yaml.load(f) + snakemake.input = Dict() + snakemake.output = Dict() + + for item in ["costs","energy"]: + snakemake.input[item] = snakemake.config['summary_dir'] + '/{name}/csvs/{item}.csv'.format(name=snakemake.config['run'],item=item) + snakemake.output[item] = snakemake.config['summary_dir'] + '/{name}/graphs/{item}.pdf'.format(name=snakemake.config['run'],item=item) + + plot_costs() + + plot_energy() diff --git a/scripts/solve_network.py b/scripts/solve_network.py new file mode 100644 index 00000000..a0f06e34 --- /dev/null +++ b/scripts/solve_network.py @@ -0,0 +1,375 @@ + +import sys + +sys.path = ["/home/vres/data/tom/lib/pypsa"] + sys.path + +import numpy as np +import pandas as pd +import logging +logger = logging.getLogger(__name__) +import gc +import os + +import pypsa + +from pypsa.descriptors import free_output_series_dataframes + +# Suppress logging of the slack bus choices +pypsa.pf.logger.setLevel(logging.WARNING) + +from vresutils.benchmark import memory_logger + + + + + +#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["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["p2"] = ["series","MW",0.,"2nd bus output","Output"] +override_component_attrs["Link"].loc["p3"] = ["series","MW",0.,"3rd bus output","Output"] + + + +def patch_pyomo_tmpdir(tmpdir): + # PYOMO should write its lp files into tmp here + import os + if not os.path.isdir(tmpdir): + os.mkdir(tmpdir) + from pyutilib.services import TempfileManager + TempfileManager.tempdir = tmpdir + +def prepare_network(n, solve_opts=None): + if solve_opts is None: + solve_opts = snakemake.config['solving']['options'] + + if 'clip_p_max_pu' in solve_opts: + for df in (n.generators_t.p_max_pu, n.storage_units_t.inflow): + df.where(df>solve_opts['clip_p_max_pu'], other=0., inplace=True) + + if solve_opts.get('load_shedding'): + n.add("Carrier", "Load") + n.madd("Generator", n.buses.index, " load", + bus=n.buses.index, + carrier='load', + sign=1e-3, # Adjust sign to measure p and p_nom in kW instead of MW + marginal_cost=1e2, # Eur/kWh + # intersect between macroeconomic and surveybased + # willingness to pay + # http://journal.frontiersin.org/article/10.3389/fenrg.2015.00055/full + p_nom=1e9 # kW + ) + + if solve_opts.get('noisy_costs'): + for t in n.iterate_components(): + #if 'capital_cost' in t.df: + # t.df['capital_cost'] += 1e1 + 2.*(np.random.random(len(t.df)) - 0.5) + if 'marginal_cost' in t.df: + t.df['marginal_cost'] += 1e-2 + 2e-3*(np.random.random(len(t.df)) - 0.5) + + for t in n.iterate_components(['Line', 'Link']): + t.df['capital_cost'] += (1e-1 + 2e-2*(np.random.random(len(t.df)) - 0.5)) * t.df['length'] + + if solve_opts.get('nhours'): + nhours = solve_opts['nhours'] + n.set_snapshots(n.snapshots[:nhours]) + n.snapshot_weightings[:] = 8760./nhours + + return n + +def add_opts_constraints(n, opts=None): + if opts is None: + opts = snakemake.wildcards.opts.split('-') + + if 'BAU' in opts: + mincaps = snakemake.config['electricity']['BAU_mincapacities'] + def bau_mincapacities_rule(model, carrier): + gens = n.generators.index[n.generators.p_nom_extendable & (n.generators.carrier == carrier)] + return sum(model.generator_p_nom[gen] for gen in gens) >= mincaps[carrier] + n.model.bau_mincapacities = pypsa.opt.Constraint(list(mincaps), rule=bau_mincapacities_rule) + + if 'SAFE' in opts: + peakdemand = (1. + snakemake.config['electricity']['SAFE_reservemargin']) * n.loads_t.p_set.sum(axis=1).max() + conv_techs = snakemake.config['plotting']['conv_techs'] + exist_conv_caps = n.generators.loc[n.generators.carrier.isin(conv_techs) & ~n.generators.p_nom_extendable, 'p_nom'].sum() + ext_gens_i = n.generators.index[n.generators.carrier.isin(conv_techs) & n.generators.p_nom_extendable] + n.model.safe_peakdemand = pypsa.opt.Constraint(expr=sum(n.model.generator_p_nom[gen] for gen in ext_gens_i) >= peakdemand - exist_conv_caps) + +def add_lv_constraint(n): + line_volume = getattr(n, 'line_volume_limit', None) + if line_volume is not None and not np.isinf(line_volume): + n.model.line_volume_constraint = pypsa.opt.Constraint( + expr=((sum(n.model.passive_branch_s_nom["Line",line]*n.lines.at[line,"length"] + for line in n.lines.index[n.lines.s_nom_extendable]) + + sum(n.model.link_p_nom[link]*n.links.at[link,"length"] + for link in n.links.index[(n.links.carrier=='DC') & + n.links.p_nom_extendable])) + <= line_volume) + ) + +def add_eps_storage_constraint(n): + if not hasattr(n, 'epsilon'): + n.epsilon = 1e-5 + fix_sus_i = n.storage_units.index[~ n.storage_units.p_nom_extendable] + n.model.objective.expr += sum(n.epsilon * n.model.state_of_charge[su, n.snapshots[0]] for su in fix_sus_i) + +def add_battery_constraints(network): + + nodes = list(network.buses.index[network.buses.carrier == "battery"]) + + def battery(model, node): + return model.link_p_nom[node + " charger"] == model.link_p_nom[node + " discharger"]*network.links.at[node + " charger","efficiency"] + + network.model.battery = pypsa.opt.Constraint(nodes, rule=battery) + + + +def add_chp_constraints(network): + + options = snakemake.config["sector"] + + if hasattr(network.links.index,"str") and network.links.index.str.contains("CHP").any(): + + #also catches central heat buses for district heating + nodes = list(network.links.index[network.links.index.str.contains("CHP electric")].str[:-len(" CHP electric")]) + + def chp_nom(model,node): + return network.links.at[node + " CHP electric","efficiency"]*options['chp_parameters']['p_nom_ratio']*model.link_p_nom[node + " CHP electric"] == network.links.at[node + " CHP heat","efficiency"]*options['chp_parameters']['p_nom_ratio']*model.link_p_nom[node + " CHP heat"] + + + network.model.chp_nom = pypsa.opt.Constraint(nodes,rule=chp_nom) + + + def backpressure(model,node,snapshot): + return options['chp_parameters']['c_m']*network.links.at[node + " CHP heat","efficiency"]*model.link_p[node + " CHP heat",snapshot] <= network.links.at[node + " CHP electric","efficiency"]*model.link_p[node + " CHP electric",snapshot] + + network.model.backpressure = pypsa.opt.Constraint(nodes,list(network.snapshots),rule=backpressure) + + + def top_iso_fuel_line(model,node,snapshot): + return model.link_p[node + " CHP heat",snapshot] + model.link_p[node + " CHP electric",snapshot] <= model.link_p_nom[node + " CHP electric"] + + network.model.top_iso_fuel_line = pypsa.opt.Constraint(nodes,list(network.snapshots),rule=top_iso_fuel_line) + + + + +def fix_branches(n, lines_s_nom=None, links_p_nom=None): + if lines_s_nom is not None and len(lines_s_nom) > 0: + for l, s_nom in lines_s_nom.iteritems(): + n.model.passive_branch_s_nom["Line", l].fix(s_nom) + if isinstance(n.opt, pypsa.opf.PersistentSolver): + n.opt.update_var(n.model.passive_branch_s_nom) + + if links_p_nom is not None and len(links_p_nom) > 0: + for l, p_nom in links_p_nom.iteritems(): + n.model.link_p_nom[l].fix(p_nom) + if isinstance(n.opt, pypsa.opf.PersistentSolver): + n.opt.update_var(n.model.link_p_nom) + +def solve_network(n, config=None, solver_log=None, opts=None): + if config is None: + config = snakemake.config['solving'] + solve_opts = config['options'] + + solver_options = config['solver'].copy() + if solver_log is None: + solver_log = snakemake.log.solver + solver_name = solver_options.pop('name') + + def run_lopf(n, allow_warning_status=False, fix_zero_lines=False, fix_ext_lines=False): + free_output_series_dataframes(n) + + if not hasattr(n, 'opt') or not isinstance(n.opt, pypsa.opf.PersistentSolver): + pypsa.opf.network_lopf_build_model(n, formulation=solve_opts['formulation']) + add_opts_constraints(n, opts) + add_lv_constraint(n) + #add_eps_storage_constraint(n) + add_battery_constraints(n) + add_chp_constraints(n) + + pypsa.opf.network_lopf_prepare_solver(n, solver_name=solver_name) + + if fix_zero_lines: + fix_lines_b = (n.lines.s_nom_opt == 0.) & n.lines.s_nom_extendable + fix_links_b = (n.links.carrier=='DC') & (n.links.p_nom_opt == 0.) & n.links.p_nom_extendable + fix_branches(n, + lines_s_nom=pd.Series(0., n.lines.index[fix_lines_b]), + links_p_nom=pd.Series(0., n.links.index[fix_links_b])) + + if fix_ext_lines: + fix_branches(n, + lines_s_nom=n.lines.loc[n.lines.s_nom_extendable, 's_nom_opt'], + links_p_nom=n.links.loc[(n.links.carrier=='DC') & n.links.p_nom_extendable, 'p_nom_opt']) + + + if not fix_ext_lines and hasattr(n.model, 'line_volume_constraint'): + + def extra_postprocessing(n, snapshots, duals): + index = list(n.model.line_volume_constraint.keys()) + cdata = pd.Series(list(n.model.line_volume_constraint.values()), + index=index) + n.line_volume_limit_dual = -cdata.map(duals).sum() + print("line volume limit dual:",n.line_volume_limit_dual) + + else: + extra_postprocessing = None + + # Firing up solve will increase memory consumption tremendously, so + # make sure we freed everything we can + gc.collect() + + #from pyomo.opt import ProblemFormat + #print("Saving model to MPS") + #n.model.write('/home/ka/ka_iai/ka_kc5996/projects/pypsa-eur/128-B-I.mps', format=ProblemFormat.mps) + #print("Model is saved to MPS") + #sys.exit() + + + status, termination_condition = \ + pypsa.opf.network_lopf_solve(n, + solver_logfile=solver_log, + solver_options=solver_options, + formulation=solve_opts['formulation'], + extra_postprocessing=extra_postprocessing + #free_memory={'pypsa'} + ) + + assert status == "ok" or allow_warning_status and status == 'warning', \ + ("network_lopf did abort with status={} " + "and termination_condition={}" + .format(status, termination_condition)) + + return status, termination_condition + + lines_ext_b = n.lines.s_nom_extendable + if lines_ext_b.any(): + # puh: ok, we need to iterate, since there is a relation + # between s/p_nom and r, x for branches. + msq_threshold = 0.01 + lines = pd.DataFrame(n.lines[['r', 'x', 'type', 'num_parallel']]) + + lines['s_nom'] = ( + np.sqrt(3) * n.lines['type'].map(n.line_types.i_nom) * + n.lines.bus0.map(n.buses.v_nom) + ).where(n.lines.type != '', n.lines['s_nom']) + + lines_ext_typed_b = (n.lines.type != '') & lines_ext_b + lines_ext_untyped_b = (n.lines.type == '') & lines_ext_b + + def update_line_parameters(n, zero_lines_below=10, fix_zero_lines=False): + if zero_lines_below > 0: + n.lines.loc[n.lines.s_nom_opt < zero_lines_below, 's_nom_opt'] = 0. + n.links.loc[(n.links.carrier=='DC') & (n.links.p_nom_opt < zero_lines_below), 'p_nom_opt'] = 0. + + if lines_ext_untyped_b.any(): + for attr in ('r', 'x'): + n.lines.loc[lines_ext_untyped_b, attr] = ( + lines[attr].multiply(lines['s_nom']/n.lines['s_nom_opt']) + ) + + if lines_ext_typed_b.any(): + n.lines.loc[lines_ext_typed_b, 'num_parallel'] = ( + n.lines['s_nom_opt']/lines['s_nom'] + ) + logger.debug("lines.num_parallel={}".format(n.lines.loc[lines_ext_typed_b, 'num_parallel'])) + + if isinstance(n.opt, pypsa.opf.PersistentSolver): + n.calculate_dependent_values() + + assert solve_opts['formulation'] == 'kirchhoff', \ + "Updating persistent solvers has only been implemented for the kirchhoff formulation for now" + + n.opt.remove_constraint(n.model.cycle_constraints) + del n.model.cycle_constraints_index + del n.model.cycle_constraints_index_0 + del n.model.cycle_constraints_index_1 + del n.model.cycle_constraints + + pypsa.opf.define_passive_branch_flows_with_kirchhoff(n, n.snapshots, skip_vars=True) + n.opt.add_constraint(n.model.cycle_constraints) + + iteration = 1 + + lines['s_nom_opt'] = lines['s_nom'] * n.lines['num_parallel'].where(n.lines.type != '', 1.) + status, termination_condition = run_lopf(n, allow_warning_status=True) + + def msq_diff(n): + lines_err = np.sqrt(((n.lines['s_nom_opt'] - lines['s_nom_opt'])**2).mean())/lines['s_nom_opt'].mean() + logger.info("Mean square difference after iteration {} is {}".format(iteration, lines_err)) + return lines_err + + min_iterations = solve_opts.get('min_iterations', 2) + max_iterations = solve_opts.get('max_iterations', 999) + + while msq_diff(n) > msq_threshold or iteration < min_iterations: + if iteration >= max_iterations: + logger.info("Iteration {} beyond max_iterations {}. Stopping ...".format(iteration, max_iterations)) + break + + update_line_parameters(n) + lines['s_nom_opt'] = n.lines['s_nom_opt'] + iteration += 1 + + status, termination_condition = run_lopf(n, allow_warning_status=True) + + update_line_parameters(n, zero_lines_below=100) + + logger.info("Starting last run with fixed extendable lines") + + # Not really needed, could also be taken out + # if 'snakemake' in globals(): + # fn = os.path.basename(snakemake.output[0]) + # n.export_to_netcdf('/home/vres/data/jonas/playground/pypsa-eur/' + fn) + + status, termination_condition = run_lopf(n, fix_ext_lines=True) + + # Drop zero lines from network + # zero_lines_i = n.lines.index[(n.lines.s_nom_opt == 0.) & n.lines.s_nom_extendable] + # if len(zero_lines_i): + # n.mremove("Line", zero_lines_i) + # zero_links_i = n.links.index[(n.links.p_nom_opt == 0.) & n.links.p_nom_extendable] + # if len(zero_links_i): + # n.mremove("Link", zero_links_i) + + + return n + +if __name__ == "__main__": + # Detect running outside of snakemake and mock snakemake for testing + if 'snakemake' not in globals(): + from vresutils.snakemake import MockSnakemake, Dict + snakemake = MockSnakemake( + wildcards=dict(network='elec', simpl='', clusters='45', lv='1.25', opts='Co2L-3H-T-H'), + input=["networks/{network}_s{simpl}_{clusters}_lv{lv}_{opts}.nc"], + output=["results/networks/s{simpl}_{clusters}_lv{lv}_{opts}-test.nc"], + log=dict(gurobi="logs/{network}_s{simpl}_{clusters}_lv{lv}_{opts}_gurobi-test.log", + python="logs/{network}_s{simpl}_{clusters}_lv{lv}_{opts}_python-test.log") + ) + + + tmpdir = snakemake.config['solving'].get('tmpdir') + if tmpdir is not None: + patch_pyomo_tmpdir(tmpdir) + + logging.basicConfig(filename=snakemake.log.python, + level=snakemake.config['logging_level']) + + with memory_logger(filename=getattr(snakemake.log, 'memory', None), interval=30.) as mem: + n = pypsa.Network(snakemake.input[0], + override_component_attrs=override_component_attrs) + + n = prepare_network(n) + n = solve_network(n) + + n.export_to_netcdf(snakemake.output[0]) + + logger.info("Maximum memory usage: {}".format(mem.mem_usage))