564 differences between costs in optimization tab and costs tab

.gitignore

@@ -167,3 +167,6 @@ cython_debug/
 
 # cached optimization results:
 ptxboa/cache
+
+# tests output:
+tests/out

ptxboa/api_calc.py

@@ -33,6 +33,7 @@ def calculate(data: CalculateDataType) -> pd.DataFrame:
             main_output_value_before_transport *= step_data["EFF"]
 
         # accumulate needed electric input
+        step_before_transport = True
         sum_el = main_output_value
         results = []
 
@@ -50,15 +51,18 @@ def calculate(data: CalculateDataType) -> pd.DataFrame:
             }
             result_process_type = df_processes.at[process_code, "result_process_type"]
 
-            main_input_value = main_output_value
-
             eff = step_data["EFF"]
-            main_output_value = main_input_value * eff
+
+            # storage efficiency must not affect main chain scaling factors:
+            if process_code not in ["EL-STR", "H2-STR"]:
+                main_input_value = main_output_value
+                main_output_value = main_input_value * eff
+
             opex_o = step_data["OPEX-O"]
 
             if not is_transport:
                 flh = step_data["FLH"]
-                liefetime = step_data["LIFETIME"]
+                lifetime = step_data["LIFETIME"]
                 capex_rel = step_data["CAPEX"]
                 opex_f = step_data["OPEX-F"]
 
@@ -72,13 +76,14 @@ def calculate(data: CalculateDataType) -> pd.DataFrame:
                     capacity = main_output_value / flh
 
                 capex = capacity * capex_rel
-                capex_ann = annuity(wacc, liefetime, capex)
+                capex_ann = annuity(wacc, lifetime, capex)
                 opex = opex_f * capacity + opex_o * main_output_value
 
                 results.append((result_process_type, process_code, "CAPEX", capex_ann))
                 results.append((result_process_type, process_code, "OPEX", opex))
 
             else:
+                step_before_transport = False
                 opex_t = step_data["OPEX-T"]
                 dist_transport = step_data["DIST"]
                 opex_ot = opex_t * dist_transport
@@ -99,14 +104,14 @@ def calculate(data: CalculateDataType) -> pd.DataFrame:
                     ]
 
                     # no FLH
-                    liefetime = sec_process_data["LIFETIME"]
+                    lifetime = sec_process_data["LIFETIME"]
                     capex = sec_process_data["CAPEX"]
                     opex_f = sec_process_data["OPEX-F"]
                     opex_o = sec_process_data["OPEX-O"]
 
                     capacity = flow_value  # no FLH
                     capex = capacity * capex
-                    capex_ann = annuity(wacc, liefetime, capex)
+                    capex_ann = annuity(wacc, lifetime, capex)
                     opex = opex_f * capacity + opex_o * flow_value
 
                     results.append(
@@ -118,9 +123,13 @@ def calculate(data: CalculateDataType) -> pd.DataFrame:
 
                     for sec_flow_code, sec_conv in sec_process_data["CONV"].items():
                         sec_flow_value = flow_value * sec_conv
-                        if sec_flow_code == "EL":
+
+                        # electricity before transport will be handled by RES step
+                        # after transport: market
+                        if sec_flow_code == "EL" and step_before_transport:
                             sum_el += sec_flow_value
-                            # TODO: in this case: no cost?
+                            # do not add SPECCOST below
+                            continue
 
                         sec_speccost = parameters["SPECCOST"][sec_flow_code]
                         sec_flow_cost = sec_flow_value * sec_speccost
@@ -142,12 +151,16 @@ def calculate(data: CalculateDataType) -> pd.DataFrame:
                 else:
                     # use market
                     speccost = parameters["SPECCOST"][flow_code]
-                    if flow_code == "EL":
+
+                    # electricity before transport will be handled by RES step
+                    # after transport: market
+                    if flow_code == "EL" and step_before_transport:
                         sum_el += flow_value
-                        # TODO: in this case: no cost?
+                        # do not add SPECCOST below
+                        continue
+
                     flow_cost = flow_value * speccost
 
-                    # TODO: not nice
                     if is_transport:
                         flow_cost = flow_cost * dist_transport
 
@@ -168,9 +181,16 @@ def calculate(data: CalculateDataType) -> pd.DataFrame:
         results = results.groupby(dim_columns).sum().reset_index()
 
         # normalization:
-        # scale so that we star twith 1 EL input,
+        # scale so that we start with 1 EL input,
         # rescale so that we have 1 unit output
-        norm_factor = sum_el / main_output_value
+        norm_factor = 1 / main_output_value
         results["values"] = results["values"] * norm_factor
 
+        # rescale again ONLY RES to account for additionally needed electricity
+        # sum_el is larger than 1.0
+        norm_factor_el = sum_el
+        idx = results["process_type"] == "Electricity generation"
+        assert idx.any()  # must have at least one entry
+        results.loc[idx, "values"] = results.loc[idx, "values"] * norm_factor_el
+
         return results

tests/test_opt.py

@@ -2,9 +2,11 @@
 """Test flh optimization."""
 
 import logging
-from json import load
+import os
+from json import dump, load
 from pathlib import Path
 from tempfile import TemporaryDirectory
+from typing import Tuple
 
 import pandas as pd
 import pypsa
@@ -15,6 +17,7 @@
 from flh_opt.api_opt import get_profiles_and_weights, optimize
 from ptxboa import DEFAULT_CACHE_DIR
 from ptxboa.api import DataHandler, PtxboaAPI
+from ptxboa.utils import annuity
 
 logging.basicConfig(level=logging.INFO)
 
@@ -217,12 +220,12 @@ def network(api) -> pypsa.Network:
 
 
 @pytest.fixture
-def network_green_iron(api) -> pypsa.Network:
+def network_green_iron(api) -> Tuple[pypsa.Network, dict, dict]:
     settings = {
         "region": "Morocco",
         "country": "Germany",
         "chain": "Green Iron (AEL)",
-        "res_gen": "Wind-PV-Hybrid",
+        "res_gen": "PV tilted",
         "scenario": "2040 (medium)",
         "secproc_co2": "Specific costs",
         "secproc_water": "Specific costs",
@@ -233,21 +236,138 @@ def network_green_iron(api) -> pypsa.Network:
     n, metadata = api.get_flh_opt_network(**settings)
     assert metadata["model_status"] == ["ok", "optimal"], "Model status not optimal"
 
-    return n, metadata
+    return n, metadata, settings
+
+
+def test_issue_564(network_green_iron, api):
+    # calculate costs from optimization tab:
+    n, metadata, settings = network_green_iron
+    res_opt = calc_aggregate_statistics(n, include_debugging_output=True)
+
+    # get costs from costs tab:
+    df_res_costs, _ = api.calculate(**settings)
+    res_costs_agg = df_res_costs.pivot_table(
+        index="process_type", columns="cost_type", values="values", aggfunc=sum
+    ).fillna(0)
+
+    res_costs_agg["total"] = res_costs_agg.sum(axis=1)
+    res_costs_agg.loc["Total"] = res_costs_agg.sum(axis=0)
+
+    # combine both sources to single df:
+    res_costs_agg.at["Electricity generation", "total_opt"] = res_opt.at[
+        "PV tilted", "Cost (USD/MWh)"
+    ]
+
+    res_costs_agg.at["Derivative production", "total_opt"] = res_opt.at[
+        "Derivative production", "Cost (USD/MWh)"
+    ]
+
+    res_costs_agg.at["Electricity and H2 storage", "total_opt"] = (
+        res_opt.at["H2 storage", "Cost (USD/MWh)"]
+        + res_opt.at["Electricity storage", "Cost (USD/MWh)"]
+    )
+
+    res_costs_agg.at["Electrolysis", "total_opt"] = res_opt.at[
+        "Electrolyzer", "Cost (USD/MWh)"
+    ]
+
+    res_costs_agg.at["Water", "total_opt"] = res_opt.at[
+        "Water supply", "Cost (USD/MWh)"
+    ]
+
+    res_costs_agg["diff"] = (
+        res_costs_agg["total_opt"] - res_costs_agg["total"]
+    ).fillna(0)
+
+    # call optimize function directly:
+    res_optimize = optimize(metadata["opt_input_data"])[0]
+
+    # write costs data to excel, and metadata to json:
+    if not os.path.exists("tests/out"):
+        os.makedirs("tests/out")
+    res_costs_agg.to_excel("tests/out/test_issue_564_res_costs_agg.xlsx")
+    res_opt.to_excel("tests/out/test_issue_564_res_opt.xlsx")
+    with open("tests/out/issue_564_metadata_optimize_input.json", "w") as f:
+        dump(metadata, f)
+    with open("tests/out/issue_564_metadata_optimize_output.json", "w") as f:
+        dump(res_optimize, f)
+
+    # extract DRI input data:
+    input_data = api.get_input_data(scenario=settings["scenario"])
+    input_data_dri = input_data.loc[
+        input_data["process_code"] == "Green iron reduction"
+    ].set_index("parameter_code")
+    wacc = input_data.loc[
+        (input_data["parameter_code"] == "WACC")
+        & (input_data["source_region_code"] == settings["region"]),
+        "value",
+    ].values[0]
+    capex = input_data_dri.at["CAPEX", "value"]
+    periods = input_data_dri.at["lifetime / amortization period", "value"]
+    opex_fix = input_data_dri.at["OPEX (fix)", "value"]
+
+    capex_ann_input = annuity(wacc, periods, capex)
+    capex_ann_opt = res_opt.at["Derivative production", "CAPEX (USD/kW)"]
+
+    # annuized capex should match:
+    assert capex_ann_input + opex_fix == pytest.approx(capex_ann_opt)
+
+    # FLH from optimization tab and optimize function output should match:
+    flh_opt_tab = res_opt.at["Derivative production", "Full load hours (h)"]
+    flh_opt_function = res_optimize["DERIV"]["FLH"] * 8760
+    assert flh_opt_tab == pytest.approx(flh_opt_function)
+
+    # for DRI, FLOW costs must be zero!
+    # because this process only comsumes electricity,
+    # and this should have zero specific cost
+    assert res_costs_agg.at["Derivative production", "FLOW"] == 0
+    assert input_data.loc["SPECCOST,,EL,,", "value"] == 0
+
+    # check power generation and consumption in optimization model:
+    power_gen = (
+        n.generators_t["p"]["PV-FIX"] * n.snapshot_weightings["generators"]
+    ).sum()
+    power_losses = (
+        -n.storage_units_t["p"]["EL_STR"] * n.snapshot_weightings["stores"]
+    ).sum()
+    power_cons_ely = (
+        n.links_t["p0"]["ELY"] * n._snapshot_weightings["generators"]
+    ).sum()
+    power_cons_deriv = (
+        n.links_t["p2"]["DERIV"] * n._snapshot_weightings["generators"]
+    ).sum()
+    power_balance = power_gen - power_losses - power_cons_ely - power_cons_deriv
+    assert power_balance == pytest.approx(0, abs=1e-6)
+
+    # check conversion coefficients of DRI:
+    # I dont write an assert statement, but the results look ok
+    df_deriv = pd.concat(
+        [n.links_t["p0"]["DERIV"], n.links_t["p1"]["DERIV"], n.links_t["p2"]["DERIV"]],
+        axis=1,
+        keys=["H2", "DRI", "EL"],
+    )
+    df_deriv["DRI_per_H2"] = df_deriv["DRI"] / df_deriv["H2"]
+    df_deriv["EL_per_DRI"] = df_deriv["EL"] / df_deriv["DRI"]
+    df_deriv = df_deriv.dropna()
+
+    # assert that differences between costs and opt tab are zero:
+    # this currently fails
+    for i in res_costs_agg["diff"]:
+        assert i == pytest.approx(0)
 
 
 def test_fix_green_iron(network_green_iron):
     """Test optimize input data: CAPEX of electricity storage should not be zero.
 
     See https://github.com/agoenergy/ptx-boa/issues/554
     """
-    n, metadata = network_green_iron
+    n, metadata, settings = network_green_iron
     assert metadata["opt_input_data"]["EL_STR"]["CAPEX_A"] != 0
 
 
 def test_output_of_final_product_is_8760mwh(network_green_iron):
     """Test that output of final process step is 8760MWh/a."""
-    n, metadata = network_green_iron
+    n, metadata, settings = network_green_iron
     res = calc_aggregate_statistics(n)
 
     assert res.at["Derivative production", "Output (MWh/a)"] == pytest.approx(8760)