Paper mill and synthetic aviation fuels

examples/36_paper_mill/36_paper_mill_mn.yaml

@@ -0,0 +1,6 @@
+name: H2Integrate_config
+system_summary: This reference paper mill plant is located in Minnesota and for its first pass, it contains paper mill plant
+  powered by grid. The system is designed to produce paper at a constant rate throughout the year.
+driver_config: driver_config.yaml
+technology_config: tech_config.yaml
+plant_config: plant_config.yaml

examples/36_paper_mill/Breakdown_cost_plots.py

@@ -0,0 +1,184 @@
+# # -*- coding: utf-8 -*-
+# """
+# Created on Fri May 15 07:38:06 2026
+
+# @author: mkoleva
+# """
+
+# import pandas as pd
+# import matplotlib.pyplot as plt
+
+# # Load Excel file
+# file_path = "Breakdown_costs_per_scenario.xlsx"
+# df = pd.read_excel(file_path, sheet_name="Sheet1", header=None)
+
+# # Scenario labels — edit however you prefer
+# scenarios = [
+#     "Paper + Pulp",
+#     "SAF with H2",
+#     "SAF with low-carbon H2",
+#     "Paper + Pulp\nSAF with H2",
+#     "Paper + Pulp\nSAF with low-carbon H2"
+# ]
+
+# # Extract cost component names
+# components = df.iloc[3:, 0].values
+
+# # Build scenario value arrays (sum of appropriate columns)
+# records = {}
+# records["Paper + Pulp"] = df.iloc[3:, [1, 2, 3]].astype(float).sum(axis=1).values
+# records["SAF with H2"] = df.iloc[3:, [3]].astype(float).sum(axis=1).values
+# records["SAF with low-carbon H2"] = df.iloc[3:, [4]].astype(float).sum(axis=1).values
+# records["Paper + Pulp\nSAF with H2"] = df.iloc[3:, [5, 6, 7]].astype(float).sum(axis=1).values
+# records["Paper + Pulp\nSAF with low-carbon H2"] = df.iloc[3:, [8, 9, 10]].astype(float).sum(axis=1).values
+
+# # Build DataFrame
+# plot_df = pd.DataFrame(records, index=components)
+# plot_df = plot_df[scenarios]  # order consistently
+
+# # Assign custom colors
+# colors = []
+# for comp in plot_df.index:
+#     if "CapEx" in comp:
+#         colors.append("navy")
+#     elif "OpEx" in comp:
+#         colors.append("orange")
+#     elif "Feedstock" in comp:
+#         colors.append("deepskyblue")
+#     elif "Taxes" in comp:
+#         colors.append("lightpink")
+#     elif "Finances" in comp:
+#         colors.append("yellowgreen")
+#     else:
+#         colors.append(None)   # Let matplotlib choose default
+
+# # Plotting
+# plt.figure(figsize=(10, 6))
+# bottom = [0] * len(scenarios)
+
+# for idx, comp in enumerate(plot_df.index):
+#     plt.bar(
+#         scenarios,
+#         plot_df.loc[comp],
+#         bottom=bottom,
+#         color=colors[idx],
+#         label=comp
+#     )
+#     bottom = [bottom[i] + plot_df.loc[comp][i] for i in range(len(scenarios))]
+
+# plt.xlabel("Scenario")
+# plt.ylabel("Cost ($/kg)")
+# plt.title("Cost Breakdown per Scenario")
+
+# # FORCE horizontal x-axis labels
+# plt.xticks(rotation=0, ha="center")
+
+# plt.legend()
+# plt.tight_layout()
+
+# plt.savefig("stacked_cost_breakdown_final.png", dpi=300)
+# plt.show()
+
+import pandas as pd
+import matplotlib.pyplot as plt
+import numpy as np
+import textwrap
+
+# -------------------------------------------------------------
+# LOAD EXCEL
+# -------------------------------------------------------------
+file_path = "Breakdown_costs_per_scenario.xlsx"
+df = pd.read_excel(file_path, sheet_name="Sheet1", header=None)
+
+# -------------------------------------------------------------
+# READ STRUCTURE
+# -------------------------------------------------------------
+scenario_row = df.iloc[0, 1:].tolist()
+product_row  = df.iloc[1, 1:].tolist()
+components   = df.iloc[2:, 0].astype(str).str.strip().tolist()
+values       = df.iloc[2:, 1:].astype(float)
+
+# -------------------------------------------------------------
+# CLEAN NANS
+# -------------------------------------------------------------
+valid = [i for i, s in enumerate(scenario_row) if str(s) != "nan"]
+scenario_row = [scenario_row[i] for i in valid]
+product_row  = [product_row[i] for i in valid]
+values       = values.iloc[:, valid]
+
+# -------------------------------------------------------------
+# MULTILINE SCENARIO LABELS (automatic wrapping)
+# -------------------------------------------------------------
+scenario_row_wrapped = [
+    "\n".join(textwrap.wrap(s, width=18)) for s in scenario_row
+]
+
+# -------------------------------------------------------------
+# BUILD MULTIINDEX
+# -------------------------------------------------------------
+tuples = list(zip(scenario_row_wrapped, product_row))
+df_plot = pd.DataFrame(values.values, index=components, columns=pd.MultiIndex.from_tuples(tuples))
+
+# -------------------------------------------------------------
+# FLATTENED PRODUCT LABELS
+# -------------------------------------------------------------
+flat_products = product_row  
+
+# -------------------------------------------------------------
+# GROUP POSITIONS FOR SCENARIO LABELS
+# -------------------------------------------------------------
+scenario_groups = {}
+for idx, scen in enumerate(scenario_row_wrapped):
+    scenario_groups.setdefault(scen, []).append(idx)
+
+x = np.arange(len(flat_products))
+
+# -------------------------------------------------------------
+# COLOR MAP
+# -------------------------------------------------------------
+color_map = {
+    "CapEx ($/kg)": "navy",
+    "OpEx ($/kg)": "orange",
+    "Feedstock ($/kg)": "deepskyblue",
+    "Taxes ($/kg)": "lightpink",
+    "Finances ($/kg)": "yellowgreen"
+}
+
+# -------------------------------------------------------------
+# PLOT
+# -------------------------------------------------------------
+plt.figure(figsize=(18, 7))
+
+bottom = np.zeros(len(x))
+
+for comp in components:
+    y = df_plot.loc[comp].values
+    plt.bar(x, y, bottom=bottom, color=color_map[comp], label=comp)
+    bottom += y
+
+# -------------------------------------------------------------
+# X‑AXIS LABELS (PRODUCT LEVEL)
+# -------------------------------------------------------------
+plt.xticks(x, flat_products, rotation=0, ha="center")
+
+# -------------------------------------------------------------
+# Y‑AXIS LABEL
+# -------------------------------------------------------------
+plt.ylabel("Levelized cost ($/kg)")
+
+plt.title("Cost Breakdown by Product and Scenario")
+
+# -------------------------------------------------------------
+# SCENARIO LABELS (CENTERED ABOVE GROUPS)
+# -------------------------------------------------------------
+ymin, ymax = plt.ylim()
+for scen, idxs in scenario_groups.items():
+    center = np.mean(idxs)
+    plt.text(center, ymax + ymax*0.04, scen,
+             ha="center", va="bottom", fontsize=11, fontweight="bold")
+
+plt.ylim(ymin, ymax * 1.25)
+
+plt.legend(title="Cost Component", bbox_to_anchor=(1.02, 1), loc="upper left")
+plt.tight_layout()
+plt.show()

examples/36_paper_mill/driver_config.yaml

@@ -0,0 +1,4 @@
+name: driver_config
+description: This analysis runs a paper mill plant and matches other examples in H2Integrate
+general:
+  folder_output: outputs

examples/36_paper_mill/paper_mill_profast_finance.py

@@ -0,0 +1,198 @@
+import pytest
+import openmdao.api as om
+from pytest import fixture
+
+from h2integrate.finances.profast_lco import ProFastLCO
+
+
+@fixture
+def profast_inputs_no1():
+    params = {
+        "analysis_start_year": 2030, #changed
+        "installation_time": 24, #24 months? 
+        "inflation_rate": 0.0, 
+        "discount_rate": 0.0948,
+        "debt_equity_ratio": 1.72,
+        "property_tax_and_insurance": 0.015, #should we use this value?
+        "total_income_tax_rate": 0.2574, #should we use this value?
+        "capital_gains_tax_rate": 0.15, #should we use this value?
+        "sales_tax_rate": 0.00,
+        "debt_interest_rate": 0.046, #should we use this value?
+        "debt_type": "Revolving debt",
+        "loan_period_if_used": 0,
+        "cash_onhand_months": 1,
+        "admin_expense": 0.00,
+    }
+    cap_items = {"depr_type": "MACRS", "depr_period": 7, "refurb": [0.0]}
+    model_inputs = {"params": params, "capital_items": cap_items}
+
+    return model_inputs
+
+
+@fixture
+def fake_filtered_tech_config():
+    tech_config = {
+        "wind": {"model_inputs": {}},
+        "solar": {"model_inputs": {}},
+        "battery": {"model_inputs": {}},
+        "natural_gas": {"model_inputs": {}},
+    }
+    return tech_config
+
+
+@fixture
+def fake_cost_dict():
+    fake_costs = {
+        "capex_adjusted_wind": 0,
+        "opex_adjusted_wind": 0,
+        "varopex_adjusted_wind": [0.0] * 0,
+        "capex_adjusted_solar": 0,
+        "opex_adjusted_solar": 0,
+        "varopex_adjusted_solar": [0.0] * 0,
+        "capex_adjusted_battery": 0,
+        "opex_adjusted_battery": 0,
+        "varopex_adjusted_battery": [0.0] * 00,
+        "capex_adjusted_natural_gas": 0,
+        "opex_adjusted_natural_gas": 0,
+        "varopex_adjusted_natural_gas": [0] * 00,
+    }
+    return fake_costs
+
+
+@pytest.mark.regression
+def test_profast_comp(profast_inputs_no1, fake_filtered_tech_config, fake_cost_dict, subtests):
+    mean_hourly_production = 34246.6 # ton/hr
+    prob = om.Problem()
+    plant_config = {
+        "plant": {
+            "plant_life": 40,
+        },
+        "finance_parameters": {"model_inputs": profast_inputs_no1},
+    }
+    pf = ProFastLCO(
+        driver_config={},
+        plant_config=plant_config,
+        tech_config=fake_filtered_tech_config,
+        commodity_type="electricity",
+        description="no1",
+    )
+    ivc = om.IndepVarComp()
+
+    ivc.add_output("rated_electricity_production", mean_hourly_production, units="kW")
+    ivc.add_output("capacity_factor", [0.9] * plant_config["plant"]["plant_life"], units="unitless")
+
+    prob.model.add_subsystem("ivc", ivc, promotes=["*"])
+    prob.model.add_subsystem("pf", pf, promotes=["rated_electricity_production", "capacity_factor"])
+    prob.setup()
+    for variable, cost in fake_cost_dict.items():
+        units = "USD" if "capex" in variable else "USD/year"
+        prob.set_val(f"pf.{variable}", cost, units=units)
+
+    prob.run_model()
+
+    lcoe = prob.get_val("pf.LCOE_no1", units="USD/(MW*h)")
+    price = prob.get_val("pf.price_electricity_no1", units="USD/(MW*h)")
+
+    wacc = prob.get_val("pf.wacc_electricity_no1", units="percent")
+    crf = prob.get_val("pf.crf_electricity_no1", units="percent")
+    profit_index = prob.get_val("pf.profit_index_electricity_no1", units="unitless")
+    irr = prob.get_val("pf.irr_electricity_no1", units="percent")
+    ipp = prob.get_val("pf.investor_payback_period_electricity_no1", units="yr")
+
+    lcoe_breakdown = prob.get_val("pf.LCOE_no1_breakdown")
+
+    with subtests.test("LCOE"):
+        assert pytest.approx(lcoe[0], rel=1e-6) == 63.8181779
+
+    with subtests.test("WACC"):
+        assert pytest.approx(wacc[0], rel=1e-6) == 0.056453864
+
+    with subtests.test("CRF"):
+        assert pytest.approx(crf[0], rel=1e-6) == 0.0674704169
+
+    with subtests.test("Profit Index"):
+        assert pytest.approx(profit_index[0], rel=1e-6) == 2.12026237778
+
+    with subtests.test("IRR"):
+        assert pytest.approx(irr[0], rel=1e-6) == 0.0948
+
+    with subtests.test("Investor payback period"):
+        assert pytest.approx(ipp[0], rel=1e-6) == 8
+
+    with subtests.test("LCOE == price"):
+        assert pytest.approx(lcoe, rel=1e-6) == price
+
+    with subtests.test("LCOE breakdown total"):
+        assert pytest.approx(lcoe_breakdown["LCOE: Total ($/kWh)"] * 1e3, rel=1e-6) == lcoe
+
+
+@pytest.mark.regression
+def test_profast_comp_coproduct(
+    profast_inputs_no1, fake_filtered_tech_config, fake_cost_dict, subtests
+):
+    mean_hourly_production = 500000.0  # kW*h
+    grid_sell_price = 63.8181779 / 1e3  # USD/(kW*h)
+    wind_sold_USD = [-1 * mean_hourly_production * 8760 * grid_sell_price] * 30
+    fake_cost_dict.update({"varopex_adjusted_wind": wind_sold_USD})
+
+    prob = om.Problem()
+    plant_config = {
+        "plant": {
+            "plant_life": 30,
+        },
+        "finance_parameters": {"model_inputs": profast_inputs_no1},
+    }
+    pf = ProFastLCO(
+        driver_config={},
+        plant_config=plant_config,
+        tech_config=fake_filtered_tech_config,
+        commodity_type="electricity",
+        description="no1",
+    )
+    ivc = om.IndepVarComp()
+    ivc.add_output("rated_electricity_production", mean_hourly_production, units="kW")
+    ivc.add_output("capacity_factor", [1.0] * plant_config["plant"]["plant_life"], units="unitless")
+
+    prob.model.add_subsystem("ivc", ivc, promotes=["*"])
+    prob.model.add_subsystem("pf", pf, promotes=["rated_electricity_production", "capacity_factor"])
+    prob.setup()
+    for variable, cost in fake_cost_dict.items():
+        units = "USD" if "capex" in variable else "USD/year"
+        prob.set_val(f"pf.{variable}", cost, units=units)
+
+    prob.run_model()
+
+    lcoe = prob.get_val("pf.LCOE_no1", units="USD/(MW*h)")
+    price = prob.get_val("pf.price_electricity_no1", units="USD/(MW*h)")
+
+    wacc = prob.get_val("pf.wacc_electricity_no1", units="percent")
+    crf = prob.get_val("pf.crf_electricity_no1", units="percent")
+    profit_index = prob.get_val("pf.profit_index_electricity_no1", units="unitless")
+    irr = prob.get_val("pf.irr_electricity_no1", units="percent")
+    ipp = prob.get_val("pf.investor_payback_period_electricity_no1", units="yr")
+
+    lcoe_breakdown = prob.get_val("pf.LCOE_no1_breakdown")
+
+    with subtests.test("LCOE"):
+        assert pytest.approx(lcoe[0], abs=1e-6) == 0
+
+    with subtests.test("WACC"):
+        assert pytest.approx(wacc[0], rel=1e-6) == 0.056453864
+
+    with subtests.test("CRF"):
+        assert pytest.approx(crf[0], rel=1e-6) == 0.0674704169
+
+    with subtests.test("Profit Index"):
+        assert pytest.approx(profit_index[0], rel=1e-6) == 2.12026237778
+
+    with subtests.test("IRR"):
+        assert pytest.approx(irr[0], rel=1e-6) == 0.0948
+
+    with subtests.test("Investor payback period"):
+        assert pytest.approx(ipp[0], rel=1e-6) == 8
+
+    with subtests.test("LCOE == price"):
+        assert pytest.approx(lcoe, rel=1e-6) == price
+
+    with subtests.test("LCOE breakdown total"):
+        assert pytest.approx(lcoe_breakdown["LCOE: Total ($/kWh)"] * 1e3, rel=1e-6) == lcoe