diff --git a/config.default.yaml b/config.default.yaml index 8cf2af31..a23cdc79 100644 --- a/config.default.yaml +++ b/config.default.yaml @@ -400,22 +400,50 @@ solving: solver: name: gurobi - threads: 4 - method: 2 # barrier - crossover: 0 - BarConvTol: 1.e-6 - Seed: 123 - AggFill: 0 - PreDual: 0 - GURO_PAR_BARDENSETHRESH: 200 - #FeasibilityTol: 1.e-6 + options: gurobi-default + + solver_options: + gurobi-default: + threads: 4 + method: 2 # barrier + crossover: 0 + BarConvTol: 1.e-6 + Seed: 123 + AggFill: 0 + PreDual: 0 + GURO_PAR_BARDENSETHRESH: 200 + seed: 10 # Consistent seed for all plattforms + gurobi-numeric-focus: + name: gurobi + NumericFocus: 3 # Favour numeric stability over speed + method: 2 # barrier + crossover: 0 # do not use crossover + BarHomogeneous: 1 # Use homogeneous barrier if standard does not converge + BarConvTol: 1.e-5 + FeasibilityTol: 1.e-4 + OptimalityTol: 1.e-4 + ObjScale: -0.5 + threads: 8 + Seed: 123 + gurobi-fallback: # Use gurobi defaults + name: gurobi + crossover: 0 + method: 2 # barrier + BarHomogeneous: 1 # Use homogeneous barrier if standard does not converge + BarConvTol: 1.e-5 + FeasibilityTol: 1.e-5 + OptimalityTol: 1.e-5 + Seed: 123 + threads: 8 + cplex-default: + threads: 4 + lpmethod: 4 # barrier + solutiontype: 2 # non basic solution, ie no crossover + barrier_convergetol: 1.e-5 + feasopt_tolerance: 1.e-6 + + cbc-default: {} # Used in CI - #name: cplex - #threads: 4 - #lpmethod: 4 # barrier - #solutiontype: 2 # non basic solution, ie no crossover - #barrier_convergetol: 1.e-5 - #feasopt_tolerance: 1.e-6 mem: 30000 #memory in MB; 20 GB enough for 50+B+I+H2; 100 GB for 181+B+I+H2 diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 93442a41..7e42838e 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -1,16 +1,9 @@ """Solve network.""" import pypsa - import numpy as np -import pandas as pd - -from pypsa.linopt import get_var, linexpr, define_constraints - -from pypsa.linopf import network_lopf, ilopf from vresutils.benchmark import memory_logger - from helper import override_component_attrs, update_config_with_sector_opts import logging @@ -64,22 +57,36 @@ def _add_land_use_constraint_m(n): n.generators.p_nom_max.clip(lower=0, inplace=True) -def prepare_network(n, solve_opts=None): +def add_co2_sequestration_limit(n, limit=200): + """Add a global constraint on the amount of Mt CO2 that can be sequestered.""" + n.carriers.loc["co2 stored", "co2_absorptions"] = -1 + n.carriers.co2_absorptions = n.carriers.co2_absorptions.fillna(0) + + limit = limit * 1e6 + for o in opts: + if not "seq" in o: continue + limit = float(o[o.find("seq")+3:]) * 1e6 + break + + n.add("GlobalConstraint", 'co2_sequestration_limit', sense="<=", constant=limit, + type="primary_energy", carrier_attribute="co2_absorptions") + + +def prepare_network(n, solve_opts=None, config=None): if 'clip_p_max_pu' in solve_opts: for df in (n.generators_t.p_max_pu, n.generators_t.p_min_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'): + # intersect between macroeconomic and surveybased willingness to pay + # http://journal.frontiersin.org/article/10.3389/fenrg.2015.00055/full 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 ) @@ -103,177 +110,138 @@ def prepare_network(n, solve_opts=None): if snakemake.config['foresight'] == 'myopic': add_land_use_constraint(n) + if n.stores.carrier.eq('co2 stored').any(): + limit = config["sector"].get("co2_sequestration_potential", 200) + add_co2_sequestration_limit(n, limit=limit) + return n def add_battery_constraints(n): + """ + Add constraint ensuring that charger = discharger: + 1 * charger_size - efficiency * discharger_size = 0 + """ + discharger_bool = n.links.index.str.contains("battery discharger") + charger_bool = n.links.index.str.contains("battery charger") - chargers_b = n.links.carrier.str.contains("battery charger") - chargers = n.links.index[chargers_b & n.links.p_nom_extendable] - dischargers = chargers.str.replace("charger", "discharger") + dischargers_ext= n.links[discharger_bool].query("p_nom_extendable").index + chargers_ext= n.links[charger_bool].query("p_nom_extendable").index - if chargers.empty or ('Link', 'p_nom') not in n.variables.index: - return + eff = n.links.efficiency[dischargers_ext].values + lhs = n.model["Link-p_nom"].loc[chargers_ext] - n.model["Link-p_nom"].loc[dischargers_ext] * eff - link_p_nom = get_var(n, "Link", "p_nom") - - lhs = linexpr((1,link_p_nom[chargers]), - (-n.links.loc[dischargers, "efficiency"].values, - link_p_nom[dischargers].values)) - - define_constraints(n, lhs, "=", 0, 'Link', 'charger_ratio') + n.model.add_constraints(lhs == 0, name="Link-charger_ratio") def add_chp_constraints(n): - electric_bool = (n.links.index.str.contains("urban central") - & n.links.index.str.contains("CHP") - & n.links.index.str.contains("electric")) - heat_bool = (n.links.index.str.contains("urban central") - & n.links.index.str.contains("CHP") - & n.links.index.str.contains("heat")) + electric = (n.links.index.str.contains("urban central") + & n.links.index.str.contains("CHP") + & n.links.index.str.contains("electric")) + heat = (n.links.index.str.contains("urban central") + & n.links.index.str.contains("CHP") + & n.links.index.str.contains("heat")) - electric = n.links.index[electric_bool] - heat = n.links.index[heat_bool] + electric_ext = n.links[electric].query("p_nom_extendable").index + heat_ext = n.links[heat].query("p_nom_extendable").index - electric_ext = n.links.index[electric_bool & n.links.p_nom_extendable] - heat_ext = n.links.index[heat_bool & n.links.p_nom_extendable] + electric_fix = n.links[electric].query("~p_nom_extendable").index + heat_fix = n.links[heat].query("~p_nom_extendable").index - electric_fix = n.links.index[electric_bool & ~n.links.p_nom_extendable] - heat_fix = n.links.index[heat_bool & ~n.links.p_nom_extendable] - - link_p = get_var(n, "Link", "p") + p = n.model["Link-p"] # dimension: [time, link] + # output ratio between heat and electricity and top_iso_fuel_line for extendable if not electric_ext.empty: + p_nom = n.model["Link-p_nom"] - link_p_nom = get_var(n, "Link", "p_nom") + lhs = (p_nom.loc[electric_ext] * (n.links.p_nom_ratio * n.links.efficiency)[electric_ext].values - + p_nom.loc[heat_ext] * n.links.efficiency[heat_ext].values) + n.model.add_constraints(lhs == 0, name='chplink-fix_p_nom_ratio') - #ratio of output heat to electricity set by p_nom_ratio - lhs = linexpr((n.links.loc[electric_ext, "efficiency"] - *n.links.loc[electric_ext, "p_nom_ratio"], - link_p_nom[electric_ext]), - (-n.links.loc[heat_ext, "efficiency"].values, - link_p_nom[heat_ext].values)) + rename = {"Link-ext": "Link"} + lhs = p.loc[:, electric_ext] + p.loc[:, heat_ext] - p_nom.rename(rename).loc[electric_ext] + n.model.add_constraints(lhs <= 0, name='chplink-top_iso_fuel_line_ext') - define_constraints(n, lhs, "=", 0, 'chplink', 'fix_p_nom_ratio') - - #top_iso_fuel_line for extendable - lhs = linexpr((1,link_p[heat_ext]), - (1,link_p[electric_ext].values), - (-1,link_p_nom[electric_ext].values)) - - define_constraints(n, lhs, "<=", 0, 'chplink', 'top_iso_fuel_line_ext') + # top_iso_fuel_line for fixed if not electric_fix.empty: + lhs = p.loc[:, electric_fix] + p.loc[:, heat_fix] + rhs = n.links.p_nom[electric_fix] + n.model.add_constraints(lhs <= rhs, name='chplink-top_iso_fuel_line_fix') - #top_iso_fuel_line for fixed - lhs = linexpr((1,link_p[heat_fix]), - (1,link_p[electric_fix].values)) - - rhs = n.links.loc[electric_fix, "p_nom"].values - - define_constraints(n, lhs, "<=", rhs, 'chplink', 'top_iso_fuel_line_fix') - + # back-pressure if not electric.empty: + lhs = (p.loc[:, heat] * (n.links.efficiency[heat] * n.links.c_b[electric].values) - + p.loc[:, electric] * n.links.efficiency[electric]) + n.model.add_constraints(lhs <= rhs, name='chplink-backpressure') - #backpressure - lhs = linexpr((n.links.loc[electric, "c_b"].values - *n.links.loc[heat, "efficiency"], - link_p[heat]), - (-n.links.loc[electric, "efficiency"].values, - link_p[electric].values)) - - define_constraints(n, lhs, "<=", 0, 'chplink', 'backpressure') - -def basename(x): - return x.split("-2")[0] def add_pipe_retrofit_constraint(n): """Add constraint for retrofitting existing CH4 pipelines to H2 pipelines.""" - gas_pipes_i = n.links.query("carrier == 'gas pipeline' and p_nom_extendable").index h2_retrofitted_i = n.links.query("carrier == 'H2 pipeline retrofitted' and p_nom_extendable").index - if h2_retrofitted_i.empty or gas_pipes_i.empty: return - - link_p_nom = get_var(n, "Link", "p_nom") - - CH4_per_H2 = 1 / n.config["sector"]["H2_retrofit_capacity_per_CH4"] - fr = "H2 pipeline retrofitted" - to = "gas pipeline" - - pipe_capacity = n.links.loc[gas_pipes_i, 'p_nom'].rename(basename) - - lhs = linexpr( - (CH4_per_H2, link_p_nom.loc[h2_retrofitted_i].rename(index=lambda x: x.replace(fr, to))), - (1, link_p_nom.loc[gas_pipes_i]) - ) - - lhs.rename(basename, inplace=True) - define_constraints(n, lhs, "=", pipe_capacity, 'Link', 'pipe_retrofit') - - -def add_co2_sequestration_limit(n, sns): - - co2_stores = n.stores.loc[n.stores.carrier=='co2 stored'].index - - if co2_stores.empty or ('Store', 'e') not in n.variables.index: + if h2_retrofitted_i.empty or gas_pipes_i.empty: return - vars_final_co2_stored = get_var(n, 'Store', 'e').loc[sns[-1], co2_stores] + p_nom = n.model["Link-p_nom"] - lhs = linexpr((1, vars_final_co2_stored)).sum() + CH4_per_H2 = 1 / n.config["sector"]["H2_retrofit_capacity_per_CH4"] + lhs = p_nom.loc[gas_pipes_i] + CH4_per_H2 * p_nom.loc[h2_retrofitted_i] + rhs = n.links.p_nom[gas_pipes_i].rename_axis("Link-ext") - limit = n.config["sector"].get("co2_sequestration_potential", 200) * 1e6 - for o in opts: - if not "seq" in o: continue - limit = float(o[o.find("seq")+3:]) * 1e6 - break - - name = 'co2_sequestration_limit' - sense = "<=" - - n.add("GlobalConstraint", name, sense=sense, constant=limit, - type=np.nan, carrier_attribute=np.nan) - - define_constraints(n, lhs, sense, limit, 'GlobalConstraint', - 'mu', axes=pd.Index([name]), spec=name) + n.model.add_constraints(lhs == rhs, name='Link-pipe_retrofit') def extra_functionality(n, snapshots): add_battery_constraints(n) add_pipe_retrofit_constraint(n) - add_co2_sequestration_limit(n, snapshots) -def solve_network(n, config, opts='', **kwargs): - solver_options = config['solving']['solver'].copy() - solver_name = solver_options.pop('name') - cf_solving = config['solving']['options'] - track_iterations = cf_solving.get('track_iterations', False) - min_iterations = cf_solving.get('min_iterations', 4) - max_iterations = cf_solving.get('max_iterations', 6) - keep_shadowprices = cf_solving.get('keep_shadowprices', True) +def solve_network(n, config, opts="", **kwargs): + set_of_options = config['solving']['solver']['options'] + solver_options = config['solving']["solver_options"][set_of_options] if set_of_options else {} + solver_name = config['solving']['solver']['name'] + cf_solving = config["solving"]["options"] + track_iterations = cf_solving.get("track_iterations", False) + min_iterations = cf_solving.get("min_iterations", 4) + max_iterations = cf_solving.get("max_iterations", 6) # add to network for extra_functionality n.config = config n.opts = opts - if cf_solving.get('skip_iterations', False): - network_lopf(n, solver_name=solver_name, solver_options=solver_options, - extra_functionality=extra_functionality, - keep_shadowprices=keep_shadowprices, **kwargs) + skip_iterations = cf_solving.get("skip_iterations", False) + if not n.lines.s_nom_extendable.any(): + skip_iterations = True + logger.info("No expandable lines found. Skipping iterative solving.") + + if skip_iterations: + status, condition = n.optimize( + solver_name=solver_name, + extra_functionality=extra_functionality, + **solver_options, + **kwargs, + ) else: - ilopf(n, solver_name=solver_name, solver_options=solver_options, - track_iterations=track_iterations, - min_iterations=min_iterations, - max_iterations=max_iterations, - extra_functionality=extra_functionality, - keep_shadowprices=keep_shadowprices, - **kwargs) + status, condition = n.optimize.optimize_transmission_expansion_iteratively( + solver_name=solver_name, + track_iterations=track_iterations, + min_iterations=min_iterations, + max_iterations=max_iterations, + extra_functionality=extra_functionality, + **solver_options, + **kwargs, + ) + + if status != "ok": + logger.warning(f"Solving status '{status}' with termination condition '{condition}'") + return n + if __name__ == "__main__": if 'snakemake' not in globals(): from helper import mock_snakemake @@ -281,10 +249,10 @@ if __name__ == "__main__": 'solve_network', simpl='', opts="", - clusters="37", + clusters="5", lv=1.0, - sector_opts='168H-T-H-B-I-A-solar+p3-dist1', - planning_horizons="2030", + sector_opts='Co2L0-3H-T-H-B-I-A-solar+p3-dist1', + planning_horizons="2050", ) logging.basicConfig(filename=snakemake.log.python, @@ -305,11 +273,9 @@ if __name__ == "__main__": overrides = override_component_attrs(snakemake.input.overrides) n = pypsa.Network(snakemake.input.network, override_component_attrs=overrides) - n = prepare_network(n, solve_opts) + n = prepare_network(n, solve_opts, config=snakemake.config) - n = solve_network(n, config=snakemake.config, opts=opts, - solver_dir=tmpdir, - solver_logfile=snakemake.log.solver) + n = solve_network(n, config=snakemake.config, opts=opts, log_fn=snakemake.log.solver) if "lv_limit" in n.global_constraints.index: n.line_volume_limit = n.global_constraints.at["lv_limit", "constant"]