Include an optional total line cost constraint

elec_s_37_lc1.25_... adds a constraint on the total line cost for an extension
by a 25%, compare with elec_s_37_lv1.25_... for the line volume limit.

`ll` is an acronym for line limit.

Snakefile

@@ -3,7 +3,7 @@ configfile: "config.yaml"
 COSTS="data/costs.csv"
 
 wildcard_constraints:
-    lv="[0-9\.]+|inf|all",
+    ll="(v|c)([0-9\.]+|opt|all)", # line limit, can be volume or cost
     simpl="[a-zA-Z0-9]*|all",
     clusters="[0-9]+m?|all",
     sectors="[+a-zA-Z0-9]+",
@@ -11,12 +11,12 @@ wildcard_constraints:
 
 rule cluster_all_elec_networks:
     input:
-        expand("networks/elec_s{simpl}_{clusters}_lv{lv}_{opts}.nc",
+        expand("networks/elec_s{simpl}_{clusters}_l{ll}_{opts}.nc",
                **config['scenario'])
 
 rule solve_all_elec_networks:
     input:
-        expand("results/networks/elec_s{simpl}_{clusters}_lv{lv}_{opts}.nc",
+        expand("results/networks/elec_s{simpl}_{clusters}_l{ll}_{opts}.nc",
                **config['scenario'])
 
 if config['enable']['prepare_links_p_nom']:
@@ -194,10 +194,10 @@ rule cluster_network:
 
 rule prepare_network:
     input: 'networks/{network}_s{simpl}_{clusters}.nc', tech_costs=COSTS
-    output: 'networks/{network}_s{simpl}_{clusters}_lv{lv}_{opts}.nc'
+    output: 'networks/{network}_s{simpl}_{clusters}_l{ll}_{opts}.nc'
     threads: 1
     resources: mem=1000
-    # benchmark: "benchmarks/prepare_network/{network}_s{simpl}_{clusters}_lv{lv}_{opts}"
+    # benchmark: "benchmarks/prepare_network/{network}_s{simpl}_{clusters}_l{ll}_{opts}"
     script: "scripts/prepare_network.py"
 
 def memory(w):
@@ -214,14 +214,14 @@ def memory(w):
         # return 4890+310 * int(w.clusters)
 
 rule solve_network:
-    input: "networks/{network}_s{simpl}_{clusters}_lv{lv}_{opts}.nc"
+    input: "networks/{network}_s{simpl}_{clusters}_l{ll}_{opts}.nc"
-    output: "results/networks/{network}_s{simpl}_{clusters}_lv{lv}_{opts}.nc"
+    output: "results/networks/{network}_s{simpl}_{clusters}_l{ll}_{opts}.nc"
     shadow: "shallow"
     log:
-        solver="logs/{network}_s{simpl}_{clusters}_lv{lv}_{opts}_solver.log",
+        solver="logs/{network}_s{simpl}_{clusters}_l{ll}_{opts}_solver.log",
-        python="logs/{network}_s{simpl}_{clusters}_lv{lv}_{opts}_python.log",
+        python="logs/{network}_s{simpl}_{clusters}_l{ll}_{opts}_python.log",
-        memory="logs/{network}_s{simpl}_{clusters}_lv{lv}_{opts}_memory.log"
+        memory="logs/{network}_s{simpl}_{clusters}_l{ll}_{opts}_memory.log"
-    benchmark: "benchmarks/solve_network/{network}_s{simpl}_{clusters}_lv{lv}_{opts}"
+    benchmark: "benchmarks/solve_network/{network}_s{simpl}_{clusters}_l{ll}_{opts}"
     threads: 4
     resources: mem=memory
     # group: "solve" # with group, threads is ignored https://bitbucket.org/snakemake/snakemake/issues/971/group-job-description-does-not-contain
@@ -230,14 +230,14 @@ rule solve_network:
 rule solve_operations_network:
     input:
         unprepared="networks/{network}_s{simpl}_{clusters}.nc",
-        optimized="results/networks/{network}_s{simpl}_{clusters}_lv{lv}_{opts}.nc"
+        optimized="results/networks/{network}_s{simpl}_{clusters}_l{ll}_{opts}.nc"
-    output: "results/networks/{network}_s{simpl}_{clusters}_lv{lv}_{opts}_op.nc"
+    output: "results/networks/{network}_s{simpl}_{clusters}_l{ll}_{opts}_op.nc"
     shadow: "shallow"
     log:
-        solver="logs/solve_operations_network/{network}_s{simpl}_{clusters}_lv{lv}_{opts}_op_solver.log",
+        solver="logs/solve_operations_network/{network}_s{simpl}_{clusters}_l{ll}_{opts}_op_solver.log",
-        python="logs/solve_operations_network/{network}_s{simpl}_{clusters}_lv{lv}_{opts}_op_python.log",
+        python="logs/solve_operations_network/{network}_s{simpl}_{clusters}_l{ll}_{opts}_op_python.log",
-        memory="logs/solve_operations_network/{network}_s{simpl}_{clusters}_lv{lv}_{opts}_op_memory.log"
+        memory="logs/solve_operations_network/{network}_s{simpl}_{clusters}_l{ll}_{opts}_op_memory.log"
-    benchmark: "benchmarks/solve_operations_network/{network}_s{simpl}_{clusters}_lv{lv}_{opts}"
+    benchmark: "benchmarks/solve_operations_network/{network}_s{simpl}_{clusters}_l{ll}_{opts}"
     threads: 4
     resources: mem=(lambda w: 5000 + 372 * int(w.clusters))
     # group: "solve_operations"
@@ -245,29 +245,29 @@ rule solve_operations_network:
 
 rule plot_network:
     input:
-        network="results/networks/{network}_s{simpl}_{clusters}_lv{lv}_{opts}.nc",
+        network="results/networks/{network}_s{simpl}_{clusters}_l{ll}_{opts}.nc",
         tech_costs=COSTS
     output:
-        only_map="results/plots/{network}_s{simpl}_{clusters}_lv{lv}_{opts}_{attr}.{ext}",
+        only_map="results/plots/{network}_s{simpl}_{clusters}_l{ll}_{opts}_{attr}.{ext}",
-        ext="results/plots/{network}_s{simpl}_{clusters}_lv{lv}_{opts}_{attr}_ext.{ext}"
+        ext="results/plots/{network}_s{simpl}_{clusters}_l{ll}_{opts}_{attr}_ext.{ext}"
     script: "scripts/plot_network.py"
 
-def summary_networks(w):
+def input_make_summary(w):
     # It's mildly hacky to include the separate costs input as first entry
     return ([COSTS] +
-            expand("results/networks/{network}_s{simpl}_{clusters}_lv{lv}_{opts}.nc",
+            expand("results/networks/{network}_s{simpl}_{clusters}_l{ll}_{opts}.nc",
                    network=w.network,
                    **{k: config["scenario"][k] if getattr(w, k) == "all" else getattr(w, k)
-                      for k in ["simpl", "clusters", "lv", "opts"]}))
+                      for k in ["simpl", "clusters", "l", "opts"]}))
 
 rule make_summary:
-    input: summary_networks
+    input: input_make_summary
-    output: directory("results/summaries/{network}_s{simpl}_{clusters}_lv{lv}_{opts}_{country}")
+    output: directory("results/summaries/{network}_s{simpl}_{clusters}_l{ll}_{opts}_{country}")
     script: "scripts/make_summary.py"
 
 rule plot_summary:
-    input: "results/summaries/{network}_s{simpl}_{clusters}_lv{lv}_{opts}_{country}"
+    input: "results/summaries/{network}_s{simpl}_{clusters}_l{ll}_{opts}_{country}"
-    output: "results/plots/summary_{summary}_{network}_s{simpl}_{clusters}_lv{lv}_{opts}_{country}.{ext}"
+    output: "results/plots/summary_{summary}_{network}_s{simpl}_{clusters}_l{ll}_{opts}_{country}.{ext}"
     script: "scripts/plot_summary.py"
 
 # Local Variables:

config.yaml

@@ -6,7 +6,9 @@ summary_dir: results
 scenario:
   sectors: [E] # ,E+EV,E+BEV,E+BEV+V2G] # [ E+EV, E+BEV, E+BEV+V2G ]
   simpl: ['']
-  lv: [1.0, 1.125, 1.25, 1.5, 2.0, 3.0, 'inf']
+  #ll: ['v1.0', 'v1.09', 'v1.125', 'v1.18', 'v1.25', 'v1.35', 'v1.5', 'v1.7', 'v2.0', 'vopt'] # line limit a 'v' prefix means volume
+  ll: ['vopt', 'copt'] #['v1.0', 'v1.125', 'v1.25', 'v1.5', 'v2.0', 'vopt'] # line limit a 'v' prefix means volume
+  #ll: ['c1.0', 'v1.125', 'v1.25', 'v1.5', 'v2.0', 'vopt'] # line limit a 'v' prefix means volume
   clusters: [37, 45, 64, 90, 128, 181, 256, 362, 512] # (2**np.r_[5.5:9:.5]).astype(int)
   opts: [Co2L-3H] #, LC-FL, LC-T, Ep-T, Co2L-T]
 

scripts/prepare_network.py

@@ -37,10 +37,56 @@ def set_line_s_max_pu(n):
     s_max_pu = np.clip(0.5 + 0.2 * (n_clusters - 37) / (200 - 37), 0.5, 0.7)
     n.lines['s_max_pu'] = s_max_pu
 
+def set_line_cost_limit(n, lc, Nyears=1.):
+    links_dc_b = n.links.carrier == 'DC'
+
+    lines_s_nom = n.lines.s_nom.where(
+        n.lines.type == '',
+        np.sqrt(3) * n.lines.num_parallel *
+        n.lines.type.map(n.line_types.i_nom) *
+        n.lines.bus0.map(n.buses.v_nom)
+    )
+
+    n.lines['capital_cost_lc'] = n.lines['capital_cost']
+    n.links['capital_cost_lc'] = n.links['capital_cost']
+    total_line_cost = ((lines_s_nom * n.lines['capital_cost_lc']).sum() +
+                       n.links.loc[links_dc_b].eval('p_nom * capital_cost_lc').sum())
+
+    if lc == 'opt':
+        costs = load_costs(Nyears, snakemake.input.tech_costs,
+                           snakemake.config['costs'], snakemake.config['electricity'])
+        update_transmission_costs(n, costs, simple_hvdc_costs=False)
+    else:
+        # Either line_volume cap or cost
+        n.lines['capital_cost'] = 0.
+        n.links.loc[links_dc_b, 'capital_cost'] = 0.
+
+    if lc == 'opt' or lc > 1.0:
+        n.lines['s_nom_min'] = lines_s_nom
+        n.lines['s_nom_extendable'] = True
+
+        n.links.loc[links_dc_b, 'p_nom_min'] = n.links.loc[links_dc_b, 'p_nom']
+        n.links.loc[links_dc_b, 'p_nom_extendable'] = True
+
+        if lc != 'opt':
+            n.line_cost_limit = lc * total_line_cost
+
+    return n
+
 def set_line_volume_limit(n, lv, Nyears=1.):
     links_dc_b = n.links.carrier == 'DC'
 
-    if np.isinf(lv):
+    lines_s_nom = n.lines.s_nom.where(
+        n.lines.type == '',
+        np.sqrt(3) * n.lines.num_parallel *
+        n.lines.type.map(n.line_types.i_nom) *
+        n.lines.bus0.map(n.buses.v_nom)
+    )
+
+    total_line_volume = ((lines_s_nom * n.lines['length']).sum() +
+                         n.links.loc[links_dc_b].eval('p_nom * length').sum())
+
+    if lv == 'opt':
         costs = load_costs(Nyears, snakemake.input.tech_costs,
                            snakemake.config['costs'], snakemake.config['electricity'])
         update_transmission_costs(n, costs, simple_hvdc_costs=True)
@@ -49,22 +95,15 @@ def set_line_volume_limit(n, lv, Nyears=1.):
         n.lines['capital_cost'] = 0.
         n.links.loc[links_dc_b, 'capital_cost'] = 0.
 
-    if lv > 1.0:
+    if lv == 'opt' or lv > 1.0:
-        lines_s_nom = n.lines.s_nom.where(
-            n.lines.type == '',
-            np.sqrt(3) * n.lines.num_parallel *
-            n.lines.type.map(n.line_types.i_nom) *
-            n.lines.bus0.map(n.buses.v_nom)
-        )
-
         n.lines['s_nom_min'] = lines_s_nom
         n.lines['s_nom_extendable'] = True
 
         n.links.loc[links_dc_b, 'p_nom_min'] = n.links.loc[links_dc_b, 'p_nom']
         n.links.loc[links_dc_b, 'p_nom_extendable'] = True
 
-        n.line_volume_limit = lv * ((lines_s_nom * n.lines['length']).sum() +
+        if lv != 'opt':
-                                    n.links.loc[links_dc_b].eval('p_nom * length').sum())
+            n.line_volume_limit = lv * total_line_volume
 
     return n
 
@@ -90,9 +129,9 @@ if __name__ == "__main__":
     if 'snakemake' not in globals():
         from vresutils.snakemake import MockSnakemake
         snakemake = MockSnakemake(
-            wildcards=dict(network='elec', simpl='', clusters='37', lv='2', opts='Co2L-3H'),
+            wildcards=dict(network='elec', simpl='', clusters='37', ll='v2', opts='Co2L-3H'),
             input=['networks/{network}_s{simpl}_{clusters}.nc'],
-            output=['networks/{network}_s{simpl}_{clusters}_lv{lv}_{opts}.nc']
+            output=['networks/{network}_s{simpl}_{clusters}_l{ll}_{opts}.nc']
         )
 
     logging.basicConfig(level=snakemake.config['logging_level'])
@@ -119,6 +158,10 @@ if __name__ == "__main__":
     # if 'Ep' in opts:
     #     add_emission_prices(n)
 
-    set_line_volume_limit(n, float(snakemake.wildcards.lv), Nyears)
+    ll_type, factor = snakemake.wildcards.ll[0], float(snakemake.wildcards.ll[1:])
+    if ll_type == 'v':
+        set_line_volume_limit(n, factor, Nyears)
+    elif ll_type == 'c':
+        set_line_cost_limit(n, factor, Nyears)
 
     n.export_to_netcdf(snakemake.output[0])

scripts/solve_network.py

@@ -43,7 +43,7 @@ def prepare_network(n, solve_opts=None):
         )
 
     if solve_opts.get('noisy_costs'):
-        for t in n.iterate_components():
+        for t in n.iterate_components(n.one_port_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:
@@ -89,6 +89,18 @@ def add_lv_constraint(n):
                     <= line_volume)
         )
 
+def add_lc_constraint(n):
+    line_cost = getattr(n, 'line_cost_limit', None)
+    if line_cost is not None and not np.isinf(line_cost):
+        n.model.line_cost_constraint = pypsa.opt.Constraint(
+            expr=((sum(n.model.passive_branch_s_nom["Line",line]*n.lines.at[line,"capital_cost_lc"]
+                        for line in n.lines.index[n.lines.s_nom_extendable]) +
+                    sum(n.model.link_p_nom[link]*n.links.at[link,"capital_cost_lc"]
+                        for link in n.links.index[(n.links.carrier=='DC') &
+                                                    n.links.p_nom_extendable]))
+                    <= line_cost)
+        )
+
 def add_eps_storage_constraint(n):
     if not hasattr(n, 'epsilon'):
         n.epsilon = 1e-5
@@ -125,6 +137,7 @@ def solve_network(n, config=None, solver_log=None, opts=None):
             pypsa.opf.network_lopf_build_model(n, formulation=solve_opts['formulation'])
             add_opts_constraints(n, opts)
             add_lv_constraint(n)
+            add_lc_constraint(n)
             # add_eps_storage_constraint(n)
 
             pypsa.opf.network_lopf_prepare_solver(n, solver_name=solver_name)