Advanced Costing Features#
The purpose of this tutorial is to introduce advanced features available through the REE Costing Framework in PrOMMiS. To learn the basics of the REE Costing Framework, please refer to the Basic Costing Features tutorial.
Introduction#
This notebook demonstrates advanced features of the REE Costing Framework. The methods described here complement the basic costing framework to add complexity and enabled more detailed cost modeling.
Note:
Calling the function `build_process_costs` on a costing block more than once is not allowed, due to the search algorithm that aggregates process and stream costs. Therefore, a new costing block will be created for each use case in this notebook where needed.
Learning Objectives#
The tutorial will take users through the following:
Importing the required tools from PrOMMiS and related repositories
Importing and using cost models from WaterTAP
Writing and importing custom, user-defined cost models
Calculating net present value and tax estimation
Standalone method for predicting process costs using economy of numbers
Standalone method for estimating cost bounds for general process types
Useful Links:
Public GitHub Repository: prommis/prommis
REE Costing Module Code: prommis/prommis
Problem Statement#
This tutorial will demonstrate costing a membrane nanofiltration process separating lithium and magnesium ions from brine in order to introduce advanced costing features:
1 Import the necessary tools#
First, import the required Python, Pyomo, IDAES, and PrOMMiS packages. These will be implemented at various stages of the tutorial.
# import pytest
import pytest
# Pyomo packages
from pyomo.environ import (
ConcreteModel,
TransformationFactory,
Constraint,
Var,
Param,
units as pyunits,
value,
)
from pyomo.network import Arc
# IDAES packages
from idaes.core import FlowsheetBlock, UnitModelBlock, UnitModelCostingBlock
from idaes.core.solvers import get_solver
from idaes.models.unit_models import Feed, Product
# PrOMMiS packages
from prommis.uky.costing.ree_plant_capcost import QGESSCosting, QGESSCostingData
# WaterTAP packages
from prommis.nanofiltration.nf_brine import (
define_feed_composition,
initialize,
solve_model,
)
from watertap.core.solvers import (
get_solver as get_watertap_solver,
) # WaterTAP has its own solver
from watertap.property_models.multicomp_aq_sol_prop_pack import (
MCASParameterBlock,
)
from watertap.unit_models.nanofiltration_DSPMDE_0D import NanofiltrationDSPMDE0D
from watertap.unit_models.pressure_changer import Pump
2 Build Process Model Flowsheet#
The REE Costing Framework supports integration with costing methods imported from WaterTAP, an external package for modeling water treatment operations such as aqueous membrane processes involving REE components. Cost accounts imported from the built-in dictionary may be attached to either unit model blocks or generic UnitModelBlock objects. In contrast, WaterTAP costing methods look for the unit model class, i.e. to import costing for a WaterTAP model the flowsheet needs to contain that WaterTAP model before calling build_process_costs. The costing framework handles passing unit model references, importing costing methods, and extracting the necessary capital and operating cost components. On the flowsheet side, users only need to define their unit models and pass them to the costing method.
2.1 Build Flowsheet#
For this example, we will build a flowsheet containing a WaterTAP unit model. We’ll see that the build_process_costs method will automatically extract the relevant costing objects.
To begin, let’s create the initial flowsheet:
# Create a Concrete Model as the top level object
m = ConcreteModel()
# Add a flowsheet object to the model
m.fs = FlowsheetBlock(dynamic=False)
2.2 Add Unit Models#
We will use the NanofiltrationDSPMDE0D unit model, a time-dependent model implementing the Donnan Steric Pore Model with Dielectric Exclusion (DSPM-DE) for nanofiltration. More information on this model may be found in the NanofiltrationDSPMDE0D documentation. The specific model is the Nanofiltration Brine model, which exists in the PrOMMiS repository:
# WaterTAP has its own solver
watertap_solver = get_watertap_solver()
# Define the property model - feed composition method is imported from the NF flowsheet
default = define_feed_composition()
m.fs.properties = MCASParameterBlock(**default)
# Add the feed and product streams
m.fs.feed = Feed(property_package=m.fs.properties)
m.fs.permeate = Product(property_package=m.fs.properties)
m.fs.retentate = Product(property_package=m.fs.properties)
# Define unit models; note that the NF model needs to be named "unit" for the imported methods to work
m.fs.pump = Pump(property_package=m.fs.properties)
m.fs.unit = NanofiltrationDSPMDE0D(property_package=m.fs.properties)
# Connect the streams and blocks
m.fs.feed_to_pump = Arc(source=m.fs.feed.outlet, destination=m.fs.pump.inlet)
m.fs.pump_to_nf = Arc(source=m.fs.pump.outlet, destination=m.fs.unit.inlet)
m.fs.nf_to_permeate = Arc(source=m.fs.unit.permeate, destination=m.fs.permeate.inlet)
m.fs.nf_to_retentate = Arc(source=m.fs.unit.retentate, destination=m.fs.retentate.inlet)
TransformationFactory("network.expand_arcs").apply_to(m)
# The methods below are imported from the NF flowsheet
# Initialize
initialize(m, watertap_solver)
# Solve
solve_model(m, watertap_solver)
# Re-solve
solve_model(m, watertap_solver)
print("Model Solved")
('Liq', 'H2O') flow_mol_phase_comp scaling factor = 0.1
('Liq', 'Li_+') flow_mol_phase_comp scaling factor = 10
('Liq', 'Mg_2+') flow_mol_phase_comp scaling factor = 10
('Liq', 'Cl_-') flow_mol_phase_comp scaling factor = 1
Cl_- adjusted: fs.feed.properties[0.0].flow_mol_phase_comp['Liq',Cl_-] was adjusted from 4.843574775634025 and fixed to 0.9219738584555871. Electroneutrality satisfied for fs.feed.properties[0.0]. Balance Result = 1.1368683772161603e-13
2025-08-19 08:01:59 [WARNING] idaes.core.util.scaling: Missing scaling factor for fs.unit.area
2025-08-19 08:01:59 [INFO] idaes.init.fs.feed: Initialization Complete.
2025-08-19 08:01:59 [INFO] idaes.init.fs.pump.control_volume: Initialization Complete
2025-08-19 08:01:59 [INFO] idaes.init.fs.pump: Initialization Complete: optimal - Optimal Solution Found
2025-08-19 08:02:01 [INFO] idaes.init.fs.unit: Initialization Complete: optimal - Optimal Solution Found
2025-08-19 08:02:01 [INFO] idaes.init.fs.permeate: Initialization Complete.
2025-08-19 08:02:01 [INFO] idaes.init.fs.retentate: Initialization Complete.
2025-08-19 08:02:01 [INFO] idaes.prommis.nanofiltration.nf_brine: Optimizing with 0 DOFs
ipopt-watertap: ipopt with user variable scaling and IDAES jacobian constraint scaling
Ipopt 3.13.2: tol=1e-08
constr_viol_tol=1e-08
acceptable_constr_viol_tol=1e-08
bound_relax_factor=0.0
honor_original_bounds=no
nlp_scaling_method=user-scaling
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt
This version of Ipopt was compiled from source code available at
https://github.com/IDAES/Ipopt as part of the Institute for the Design of
Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.
This version of Ipopt was compiled using HSL, a collection of Fortran codes
for large-scale scientific computation. All technical papers, sales and
publicity material resulting from use of the HSL codes within IPOPT must
contain the following acknowledgement:
HSL, a collection of Fortran codes for large-scale scientific
computation. See http://www.hsl.rl.ac.uk.
******************************************************************************
This is Ipopt version 3.13.2, running with linear solver ma27.
Number of nonzeros in equality constraint Jacobian...: 1140
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 528
Total number of variables............................: 378
variables with only lower bounds: 218
variables with lower and upper bounds: 144
variables with only upper bounds: 0
Total number of equality constraints.................: 378
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 0.0000000e+00 1.53e+01 1.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 0.0000000e+00 5.49e-01 9.66e+02 -1.0 1.03e-02 - 9.90e-01 9.64e-01h 1
2 0.0000000e+00 6.94e-07 8.99e+01 -1.0 4.68e-04 - 1.00e+00 1.00e+00h 1
3 0.0000000e+00 1.42e-14 6.61e-04 -1.0 6.32e-07 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 3
(scaled) (unscaled)
Objective...............: 0.0000000000000000e+00 0.0000000000000000e+00
Dual infeasibility......: 0.0000000000000000e+00 0.0000000000000000e+00
Constraint violation....: 1.4210854715202004e-14 1.4210854715202004e-14
Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00
Overall NLP error.......: 1.4210854715202004e-14 1.4210854715202004e-14
Number of objective function evaluations = 4
Number of objective gradient evaluations = 4
Number of equality constraint evaluations = 4
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 4
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 3
Total CPU secs in IPOPT (w/o function evaluations) = 0.003
Total CPU secs in NLP function evaluations = 0.000
EXIT: Optimal Solution Found.
2025-08-19 08:02:01 [INFO] idaes.prommis.nanofiltration.nf_brine: Optimizing with 0 DOFs
ipopt-watertap: ipopt with user variable scaling and IDAES jacobian constraint scaling
Ipopt 3.13.2: tol=1e-08
constr_viol_tol=1e-08
acceptable_constr_viol_tol=1e-08
bound_relax_factor=0.0
honor_original_bounds=no
nlp_scaling_method=user-scaling
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt
This version of Ipopt was compiled from source code available at
https://github.com/IDAES/Ipopt as part of the Institute for the Design of
Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.
This version of Ipopt was compiled using HSL, a collection of Fortran codes
for large-scale scientific computation. All technical papers, sales and
publicity material resulting from use of the HSL codes within IPOPT must
contain the following acknowledgement:
HSL, a collection of Fortran codes for large-scale scientific
computation. See http://www.hsl.rl.ac.uk.
******************************************************************************
This is Ipopt version 3.13.2, running with linear solver ma27.
Number of nonzeros in equality constraint Jacobian...: 1140
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 528
Total number of variables............................: 378
variables with only lower bounds: 218
variables with lower and upper bounds: 144
variables with only upper bounds: 0
Total number of equality constraints.................: 378
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 0.0000000e+00 1.53e+01 1.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 0.0000000e+00 5.49e-01 9.66e+02 -1.0 1.03e-02 - 9.90e-01 9.64e-01h 1
2 0.0000000e+00 6.94e-07 8.99e+01 -1.0 4.68e-04 - 1.00e+00 1.00e+00h 1
3 0.0000000e+00 1.42e-14 6.61e-04 -1.0 6.32e-07 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 3
(scaled) (unscaled)
Objective...............: 0.0000000000000000e+00 0.0000000000000000e+00
Dual infeasibility......: 0.0000000000000000e+00 0.0000000000000000e+00
Constraint violation....: 1.4210854715202004e-14 1.4210854715202004e-14
Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00
Overall NLP error.......: 1.4210854715202004e-14 1.4210854715202004e-14
Number of objective function evaluations = 4
Number of objective gradient evaluations = 4
Number of equality constraint evaluations = 4
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 4
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 3
Total CPU secs in IPOPT (w/o function evaluations) = 0.003
Total CPU secs in NLP function evaluations = 0.000
EXIT: Optimal Solution Found.
Model Solved
Let’s explore the unit we just built to check for costing objects:
for i in m.fs.unit.component_data_objects(descend_into=True):
if "cost" in str(i):
print(i)
print("Search complete")
Search complete
As expected, there are no costing objects on the model yet. We don’t need to add this ourselves, as the REE Costing Framework will call the correct objects from WaterTAP automatically when we build the plant costs.
2.3 Build Process Costs#
Next, let’s define some process-wide parameters such as the product and resource rates and build process costs. In this case, the costing methods expects pure and mixed product rates; in this scenario, there is no saleable product since the retentate requires further processing, so we’ll use the retentate flow as a “treated stream” product and tell the costing that it can’t be sold (price = 0). There is still a waste disposal cost for the permeate and cost of chemicals per gallon (arbitrarily chosen for this example), which are included below:
# Define resources
m.fs.water = Var([0], initialize=1000, units=pyunits.gallon / pyunits.hr)
m.fs.water.fix()
m.fs.chemicals = Var([0], initialize=20, units=pyunits.gallon / pyunits.hr)
m.fs.chemicals.fix()
m.fs.waste = Var([0], initialize=20, units=pyunits.gallon / pyunits.hr)
m.fs.waste_constraint = Constraint(
expr=(
m.fs.waste[0]
== pyunits.convert(
m.fs.unit.feed_side.properties_out[0].flow_vol_phase["Liq"],
to_units=pyunits.gal / pyunits.hr,
)
)
)
# Define product rate - required for costing, but given price = 0 in this scenario (not a saleable product)
m.fs.product = Var([0], initialize=20, units=pyunits.gallon / pyunits.hr)
m.fs.product_constraint = Constraint(
expr=(
m.fs.product[0]
== pyunits.convert(
m.fs.unit.permeate_side[0, 1].flow_vol_phase["Liq"],
to_units=pyunits.gal / pyunits.hr,
)
)
)
# Attach a flowsheet costing block
m.fs.costing_1 = QGESSCosting()
# Build process costs
m.fs.costing_1.build_process_costs(
Lang_factor=2.97,
fixed_OM=True,
# treated stream is not saleable, so pass dummy inputs
pure_product_output_rates={"treated_stream": m.fs.product[0]},
mixed_product_output_rates={"treated_stream": m.fs.product[0]},
sale_prices={"treated_stream": 0 * pyunits.USD_2021 / pyunits.gal},
variable_OM=True,
resources=[
"water",
"chemicals",
"waste",
],
rates=[m.fs.water, m.fs.chemicals, m.fs.waste],
prices={
"chemicals": 1.00 * pyunits.USD_2021 / pyunits.gallon,
"waste": 0.35 * pyunits.USD_2021 / pyunits.gallon,
},
watertap_blocks=[
m.fs.unit,
],
CE_index_year="2021",
)
QGESSCostingData.costing_initialization(m.fs.costing_1)
QGESSCostingData.initialize_fixed_OM_costs(m.fs.costing_1)
QGESSCostingData.initialize_variable_OM_costs(m.fs.costing_1)
# Solve and display results
solver = get_solver()
solver.solve(m, tee=False)
# m.fs.costing_1.report() # uncomment to view cost results
assert m.fs.costing_1.total_BEC.value == pytest.approx(0.0035218, rel=1e-4)
assert m.fs.costing_1.total_installation_cost.value == pytest.approx(
0.0069380, rel=1e-4
)
assert m.fs.costing_1.total_plant_cost.value == pytest.approx(0.010460, rel=1e-4)
assert m.fs.costing_1.total_fixed_OM_cost.value == pytest.approx(5.1756, rel=1e-4)
assert m.fs.costing_1.total_sales_revenue.value == pytest.approx(0, abs=1e-8)
assert m.fs.costing_1.total_variable_OM_cost[0].value == pytest.approx(3.9968, rel=1e-4)
assert m.fs.costing_1.plant_overhead_cost[0].value == pytest.approx(1.0351, rel=1e-4)
assert m.fs.costing_1.watertap_fixed_costs.value == pytest.approx(0.000352180, rel=1e-4)
WARNING: model contains export suffix 'scaling_factor' that contains 27
component keys that are not exported as part of the NL file. Skipping.
WARNING: model contains export suffix 'scaling_factor' that contains 5 keys
that are not Var, Constraint, Objective, or the model. Skipping.
Checking some of the model objects:
for i in m.fs.costing_1.component_data_objects(Var, descend_into=True):
if "watertap" in str(i) or "BEC" in str(i):
print(i.name, "has a value of ", value(i))
print()
for i in m.fs.costing_1.component_data_objects(Constraint, descend_into=True):
if "watertap" in str(i.expr) or "BEC" in str(i.expr):
print(i.expr)
print()
fs.costing_1.total_BEC has a value of 0.0035218040126026163
fs.costing_1.watertap_fixed_costs has a value of 0.00035218040126027783
fs.costing_1.total_BEC == 1.1739346708671863e-06*(MUSD_2021/USD_2018)*fs.unit.costing.capital_cost
fs.costing_1.total_installation_cost == (fs.costing_1.Lang_factor - 1)*fs.costing_1.total_BEC
fs.costing_1.total_plant_cost == fs.costing_1.total_BEC + fs.costing_1.total_installation_cost + fs.costing_1.other_plant_costs
fs.costing_1.watertap_fixed_costs == 1.1739346708671863e-06*(MUSD_2021/USD_2018)*fs.unit.costing.fixed_operating_cost
fs.costing_1.total_fixed_OM_cost == fs.costing_1.annual_labor_cost + fs.costing_1.maintenance_and_material_cost + fs.costing_1.quality_assurance_and_control_cost + fs.costing_1.admin_and_support_labor_cost + fs.costing_1.sales_patenting_and_research_cost + fs.costing_1.property_taxes_and_insurance_cost + fs.costing_1.other_fixed_costs + fs.costing_1.watertap_fixed_costs + fs.costing_1.custom_fixed_costs
We see that the model has a variable called watertap_fixed_costs. This is a placeholder variable to catch operating costs returned by WaterTAP. There is no watertap_variable_costs variable, as no current WaterTAP models contains variable operating costs. The expressions for total_BEC and watertap_fixed_costs expressions each include terms for the NanofiltrationDSPMDE0D costs.
Multiple WaterTAP models may be passed to the argument watertap_blocks, and the REE Costing Framework will automatically extract all capital and operating cost components.
3 User-Defined Custom Costing#
Suppose we want to add costing for a custom unit model to the flowsheet modeling a vessel to wash and extract the permeate stream. For this example, we’ll assume that this makes the stream saleable. We know the volume of the vessel and water injection rate.
The UnitModelCostingBlock takes several arguments, including costing_method which points to an external class containing a function that defines the cost model. The argument costing_method_arguments aligns with the inputs of that function. Users may write their own functions, which is demonstrated below.
3.1 Writing a Custom Cost Model#
Suppose we have written our own costing model for the vessel. The equations below have been arbitrarily chosen for demonstrative purposes.
The capital cost (in dollars) assumes a six-tenths power law based on vessel volume:
Capital_Cost = Reference_Cost * (Volume / Reference_Volume) ** 0.6
Fixed operating costs (in dollars per year) are a set percentage of the capital cost:
Fixed_Costs = 0.05 * Capital_Cost / year
Variable operating costs (in dollars per year) are a function of the water injection rate:
Variable_Operating_Costs = ($.00190/gal) * Water_Injection_Rate
We’ll add arguments for the vessel volume, material, water injection rate, and number of units. We’ll also allow a string for the cost year, which will set the currency units for the cost variables.
The custom model is formulated below:
# Import necessary IDAES functions to declare a new process block class
from idaes.core import declare_process_block_class, FlowsheetCostingBlockData
# Import a Pyomo function for variable bound ranges
from pyomo.environ import NonNegativeReals
@declare_process_block_class(
"CustomCosting"
) # a new class should always open with a declaration decorator
class CustomCostingData(
FlowsheetCostingBlockData
): # the class inherits all properties of the IDAES FlowsheetCostingBlock class
# Register custom currency units used in the custom costing model
pyunits.load_definitions_from_strings(
[
"USD_custom = 500/700 * USD_CE500",
]
) # CEPCI = 700, USD_CE500 is a placeholder that cancels out in the unit conversion
# Define global parameters - base cost year and period
def build_global_params(self):
# Set the base year for all costs
self.base_currency = pyunits.USD_2022
# Set a base period for all operating costs
self.base_period = pyunits.year
# Custom costing method
def cost_custom_vessel(
blk, # when the costing block is built, blk will be the costing block itself
volume_per_unit=1000 * pyunits.m**3,
material="carbonsteel",
water_injection_rate_per_unit=1 * pyunits.m**3 / pyunits.s,
number_of_units=1,
):
# Make parameter for number of units
blk.number_of_units = Param(
initialize=number_of_units, mutable=True, units=pyunits.dimensionless
)
# Define the bare erected cost (BEC) per unit
material_factor_dict = (
{ # material factors for each valid choice of shell material
"carbonsteel": 1,
"stainlessteel": 1.5,
"aluminum": 2,
}
)
blk.capital_cost_per_unit = Var(
initialize=1000,
units=blk.costing_package.base_currency,
domain=NonNegativeReals,
bounds=(0, None),
)
@blk.Constraint()
def capital_cost_per_unit_eq(blk):
# cost equation CAPITAL_COST = REF_COST * (VOLUME / REF_VOLUME)**0.6
ref_cost = (
10000 * pyunits.USD_custom
) # reference cost is in reference currency units
ref_volume = 1000 * pyunits.m**3
return blk.capital_cost_per_unit == pyunits.convert(
material_factor_dict[material]
* ref_cost
* (volume_per_unit / ref_volume) ** 0.6,
to_units=blk.costing_package.base_currency, # convert to costing block base currency
)
# Create a variable capital_cost that the REE Costing Framework can look for
blk.capital_cost = Var(
initialize=1000,
units=blk.costing_package.base_currency, # define in costing block base currency
domain=NonNegativeReals,
bounds=(0, None),
)
@blk.Constraint()
def capital_cost_constraint(blk):
return blk.capital_cost == blk.capital_cost_per_unit * blk.number_of_units
# Define the fixed costs
# Create a variable fixed_operating_cost that the REE Costing Framework can look for
blk.fixed_operating_cost = Var(
initialize=1000,
units=blk.costing_package.base_currency
/ blk.costing_package.base_period, # define in costing block base currency
domain=NonNegativeReals,
bounds=(0, None),
)
# Set fixed opex = 5% of capex
# This is an arbitrary choice for this example, and not a general rule
# In practice users should create Param and Var as needed for their model
blk.fixed_opex_coefficient = Param(
initialize=0.05, mutable=True, units=pyunits.dimensionless
)
@blk.Constraint()
def fixed_operating_cost_constraint(blk):
return blk.fixed_operating_cost == pyunits.convert(
blk.fixed_opex_coefficient * blk.capital_cost / pyunits.year,
to_units=blk.costing_package.base_currency
/ blk.costing_package.base_period,
)
# Define the variable costs
blk.variable_operating_cost_per_unit = Var(
initialize=1000,
units=blk.costing_package.base_currency
/ blk.costing_package.base_period, # define in costing block base currency / time
domain=NonNegativeReals,
bounds=(0, None),
)
# Set variable opex = $0.0019 per gallon of water injected using the reference year (USD_custom)
# This is an arbitrary choice for this example, and not a general rule
blk.variable_opex_price = Param(
initialize=0.00190,
units=pyunits.USD_custom
/ pyunits.gal, # define in reference currency units
mutable=True,
)
@blk.Constraint()
def variable_operating_cost_per_unit_eq(blk):
return blk.variable_operating_cost_per_unit == pyunits.convert(
blk.variable_opex_price * water_injection_rate_per_unit,
to_units=blk.costing_package.base_currency
/ blk.costing_package.base_period, # define in costing block base currency / time
)
# Create a variable variable_operating_cost that the REE Costing Framework can look for
blk.variable_operating_cost = Var(
initialize=1000,
units=blk.costing_package.base_currency
/ blk.costing_package.base_period, # define in costing block base currency / time
domain=NonNegativeReals,
bounds=(0, None),
)
@blk.Constraint()
def variable_operating_cost_constraint(blk):
return (
blk.variable_operating_cost
== blk.variable_operating_cost_per_unit * blk.number_of_units
)
There are a few important details in the model above. First, the decorator @declare_process_block_class() tells the Python interpreter to register the new class as an importable object with the name “CustomCosting”. The class inherits all properties of the FlowsheetCostingBlock class by calling its corresponding data class.
The first line inside the custom class defines a new cost unit USD_custom with a Chemical Engineering Plant Cost Index (CEPCI) value of 700. USD_CE500 is a reference unit container enabling conversion between currency units. The method build_global_params defines the base cost year and period for the model. Since the class inherits all properties of the FlowsheetCostingBlockData class, the final model will return cost variables in units of USD_2021 and any time units will be year or a (annum, the two are equivalent in Pyomo unit container syntax).
Then, the method cost_custom_vessel defines all cost equations. Each of the three cost variables (capital, fixed operating, variable operating) is built by defining a cost per unit with appropriate inputs, and then a total cost based on the number of units. The naming convention is important, since the REE Costing Framework will look for variables named capital_cost, fixed_operating_cost, and variable_operating_cost when building the process cost equations.
Finally, note that the currency units are set as costing_package.base_currency everywhere in the custom model. The object costing_package is the main flowsheet costing class set by the user, so the model will automatically use those default cost units. The examples below will demonstrate this further.
3.2 Create a Costing Block#
Let’s create a new flowsheet using our custom costing model. We’ll use a UnitModelBlock for the vessel itself and attach it to the custom class object. The vessel should be large enough to hold 3 days worth of product flow, and the water injection rate will be 5% of the product flow (these are arbitrarily chosen for this example):
# Define vessel model
m.fs.custom_vessel = UnitModelBlock()
m.fs.custom_vessel.volume = Var(initialize=5000, units=pyunits.m**3)
m.fs.custom_vessel.volume_constraint = Constraint(
expr=(
m.fs.custom_vessel.volume
== pyunits.convert(m.fs.product[0] * 3 * pyunits.day, to_units=pyunits.m**3)
)
)
m.fs.custom_vessel.water_injection_rate = Var(
initialize=0.1, units=pyunits.m**3 / pyunits.s
)
m.fs.custom_vessel.water_injection_rate_constraint = Constraint(
expr=(
m.fs.custom_vessel.water_injection_rate
== pyunits.convert(m.fs.product[0] * 5 / 100, to_units=pyunits.m**3 / pyunits.s)
)
)
# Attach a flowsheet costing block
m.fs.costing_2 = QGESSCosting()
# Add costing
m.fs.custom_vessel.costing = UnitModelCostingBlock(
flowsheet_costing_block=m.fs.costing_2,
costing_method=CustomCostingData.cost_custom_vessel, # this is the custom method we wrote above
costing_method_arguments={
"volume_per_unit": m.fs.custom_vessel.volume,
"material": "carbonsteel",
"water_injection_rate_per_unit": m.fs.custom_vessel.water_injection_rate,
"number_of_units": 1,
# no CE_index_year argument, since our class doesn't accept this as an argument
# the cost variables will be returned in units of base_currency = USD_2021
},
)
# Solve and display results (uncomment line to display)
solver.solve(m, tee=True)
# m.fs.custom_vessel.display() # uncomment to view cost results
assert m.fs.custom_vessel.costing.capital_cost.value == pytest.approx(294.011, rel=1e-4)
assert m.fs.custom_vessel.costing.fixed_operating_cost.value == pytest.approx(
14.7006, rel=1e-4
)
assert m.fs.custom_vessel.costing.variable_operating_cost.value == pytest.approx(
8.49301, rel=1e-4
)
WARNING: model contains export suffix 'scaling_factor' that contains 27
component keys that are not exported as part of the NL file. Skipping.
WARNING: model contains export suffix 'scaling_factor' that contains 5 keys
that are not Var, Constraint, Objective, or the model. Skipping.
Ipopt 3.13.2: nlp_scaling_method=gradient-based
tol=1e-06
max_iter=200
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt
This version of Ipopt was compiled from source code available at
https://github.com/IDAES/Ipopt as part of the Institute for the Design of
Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.
This version of Ipopt was compiled using HSL, a collection of Fortran codes
for large-scale scientific computation. All technical papers, sales and
publicity material resulting from use of the HSL codes within IPOPT must
contain the following acknowledgement:
HSL, a collection of Fortran codes for large-scale scientific
computation. See http://www.hsl.rl.ac.uk.
******************************************************************************
This is Ipopt version 3.13.2, running with linear solver ma27.
Number of nonzeros in equality constraint Jacobian...: 1209
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 529
Total number of variables............................: 410
variables with only lower bounds: 241
variables with lower and upper bounds: 144
variables with only upper bounds: 0
Total number of equality constraints.................: 410
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 0.0000000e+00 1.60e+06 1.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 0.0000000e+00 1.59e+06 3.04e+04 -1.0 4.01e+04 - 6.94e-03 7.56e-03h 1
2 0.0000000e+00 1.59e+06 3.05e+04 -1.0 3.92e+04 - 1.44e-02 7.65e-05h 1
3r 0.0000000e+00 1.59e+06 9.99e+02 4.4 0.00e+00 - 0.00e+00 3.83e-07R 2
4r 0.0000000e+00 1.59e+06 6.75e+04 4.4 4.03e+05 - 2.95e-03 8.45e-06f 1
5r 0.0000000e+00 1.54e+06 1.19e+05 3.7 4.03e+05 - 5.11e-01 1.53e-02f 1
6r 0.0000000e+00 1.49e+06 1.16e+05 2.3 5.03e+05 - 1.66e-01 2.37e-02f 1
Error in an AMPL evaluation. Run with "halt_on_ampl_error yes" to see details.
Warning: Cutting back alpha due to evaluation error
Error in an AMPL evaluation. Run with "halt_on_ampl_error yes" to see details.
Warning: Cutting back alpha due to evaluation error
7r 0.0000000e+00 1.49e+06 1.16e+05 2.3 2.57e+06 -2.0 1.36e-03 1.36e-04f 3
8r 0.0000000e+00 1.28e+06 1.05e+05 2.3 3.31e+03 -0.7 1.39e-01 9.76e-02f 1
Error in an AMPL evaluation. Run with "halt_on_ampl_error yes" to see details.
Warning: Cutting back alpha due to evaluation error
Error in an AMPL evaluation. Run with "halt_on_ampl_error yes" to see details.
Warning: Cutting back alpha due to evaluation error
Error in an AMPL evaluation. Run with "halt_on_ampl_error yes" to see details.
Warning: Cutting back alpha due to evaluation error
Error in an AMPL evaluation. Run with "halt_on_ampl_error yes" to see details.
Warning: Cutting back alpha due to evaluation error
9r 0.0000000e+00 1.25e+06 1.03e+05 2.3 1.73e+04 -1.1 7.41e-03 1.85e-02f 5
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10r 0.0000000e+00 1.24e+06 1.02e+05 2.3 3.86e+03 -0.7 2.13e-01 5.42e-03f 1
Error in an AMPL evaluation. Run with "halt_on_ampl_error yes" to see details.
Warning: Cutting back alpha due to evaluation error
Error in an AMPL evaluation. Run with "halt_on_ampl_error yes" to see details.
Warning: Cutting back alpha due to evaluation error
Error in an AMPL evaluation. Run with "halt_on_ampl_error yes" to see details.
Warning: Cutting back alpha due to evaluation error
Error in an AMPL evaluation. Run with "halt_on_ampl_error yes" to see details.
Warning: Cutting back alpha due to evaluation error
11r 0.0000000e+00 1.24e+06 1.02e+05 2.3 7.46e+05 -1.2 2.18e-04 2.07e-04f 5
12r 0.0000000e+00 8.47e+05 7.67e+04 2.3 2.82e+03 -0.8 2.15e-01 2.46e-01f 1
13r 0.0000000e+00 4.47e+05 5.50e+04 2.3 1.11e+03 -0.3 1.37e-01 3.80e-01f 1
14r 0.0000000e+00 2.97e+05 3.67e+04 2.3 7.63e+03 -0.8 3.71e-01 2.65e-01f 1
15r 0.0000000e+00 1.90e+05 2.28e+04 2.3 2.21e+03 -0.4 4.38e-01 3.06e-01f 1
16r 0.0000000e+00 1.65e+05 1.96e+04 2.3 7.39e+03 -0.9 3.01e-01 1.17e-01f 1
17r 0.0000000e+00 6.20e+04 7.74e+03 2.3 2.49e+03 -0.4 7.20e-01 5.79e-01f 1
18r 0.0000000e+00 5.36e+04 7.30e+03 1.6 1.09e+02 -0.9 6.31e-01 1.32e-01f 1
19r 0.0000000e+00 1.90e+04 4.40e+03 1.6 2.43e+02 -1.4 7.50e-01 6.37e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20r 0.0000000e+00 4.02e+03 3.27e+03 1.6 1.81e+03 - 8.67e-01 7.87e-01f 1
21r 0.0000000e+00 2.49e+03 3.03e+04 0.9 9.33e+02 - 8.57e-01 3.81e-01f 1
22r 0.0000000e+00 1.34e+03 3.34e+04 0.9 2.86e+03 - 9.77e-01 4.62e-01f 1
23r 0.0000000e+00 1.33e+02 6.87e+03 0.9 1.34e+03 - 1.00e+00 9.49e-01f 1
24r 0.0000000e+00 1.31e+02 1.40e+03 0.9 1.07e+02 - 1.00e+00 1.00e+00f 1
25r 0.0000000e+00 1.30e+02 1.26e+01 0.9 2.42e+01 - 1.00e+00 1.00e+00h 1
26r 0.0000000e+00 5.36e+01 4.22e+03 -1.2 1.88e+03 - 7.63e-01 6.90e-01f 1
27r 0.0000000e+00 2.53e+01 7.81e+03 -1.2 2.00e+02 - 7.74e-01 5.28e-01f 1
28r 0.0000000e+00 1.30e+01 6.46e+03 -1.2 5.59e+02 - 5.89e-01 4.66e-01f 1
29r 0.0000000e+00 7.68e+00 6.48e+03 -1.2 6.69e+02 - 1.00e+00 6.05e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
30r 0.0000000e+00 1.18e+01 3.06e+03 -1.2 1.10e+03 - 8.32e-01 7.24e-01f 1
31r 0.0000000e+00 1.91e+01 5.33e+03 -1.2 7.70e+02 - 1.00e+00 1.00e+00f 1
32r 0.0000000e+00 1.81e+01 9.62e+01 -1.2 1.17e+02 - 1.00e+00 1.00e+00f 1
33r 0.0000000e+00 1.80e+01 1.62e-01 -1.2 5.18e-01 - 1.00e+00 1.00e+00h 1
34r 0.0000000e+00 2.68e+01 2.39e+03 -2.8 9.13e+02 - 5.91e-01 6.63e-01f 1
35r 0.0000000e+00 2.49e+01 1.71e+03 -2.8 1.99e+02 - 7.72e-01 6.91e-01f 1
36r 0.0000000e+00 2.52e+01 1.10e+03 -2.8 7.17e+01 - 1.00e+00 7.49e-01f 1
37r 0.0000000e+00 2.77e+01 2.76e+02 -2.8 3.70e+01 - 1.00e+00 1.00e+00f 1
38r 0.0000000e+00 2.82e+01 8.92e+00 -2.8 1.80e+00 - 1.00e+00 1.00e+00h 1
39r 0.0000000e+00 2.82e+01 7.24e-03 -2.8 1.31e-01 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
40r 0.0000000e+00 3.31e+01 1.66e+02 -4.2 8.47e+01 - 8.53e-01 8.87e-01f 1
41r 0.0000000e+00 3.33e+01 3.46e+01 -4.2 3.38e-01 -1.9 1.00e+00 9.51e-01f 1
42r 0.0000000e+00 3.38e+01 2.74e+02 -4.2 1.01e+00 -2.4 6.10e-01 1.00e+00f 1
43r 0.0000000e+00 3.38e+01 1.43e+02 -4.2 3.00e+00 -2.8 8.62e-01 5.94e-01f 1
44r 0.0000000e+00 3.37e+01 6.14e+01 -4.2 8.73e+00 -3.3 3.85e-01 1.27e-01f 1
45r 0.0000000e+00 3.24e+01 2.14e+02 -4.2 2.44e+01 -3.8 1.00e+00 5.17e-01f 1
46r 0.0000000e+00 3.06e+01 1.42e+02 -4.2 5.84e+01 -4.3 2.86e-01 3.00e-01f 1
47r 0.0000000e+00 3.02e+01 5.73e+02 -4.2 1.10e+02 -4.7 1.00e+00 3.10e-02h 1
48r 0.0000000e+00 1.80e+01 2.53e+02 -4.2 1.65e+02 -5.2 1.00e+00 7.22e-01f 1
49r 0.0000000e+00 1.70e+01 2.42e+02 -4.2 1.27e+02 -5.7 1.00e+00 1.24e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
50r 0.0000000e+00 1.28e+01 3.21e+02 -4.2 1.00e+02 -6.2 1.00e+00 5.26e-01h 1
51r 0.0000000e+00 1.25e+01 1.09e+03 -4.2 2.27e+02 -6.7 6.82e-01 7.12e-02h 1
52r 0.0000000e+00 1.04e+01 5.57e+02 -4.2 7.42e+02 -7.1 6.05e-01 5.49e-01h 1
53r 0.0000000e+00 9.50e+00 4.56e+02 -4.2 3.22e+03 -7.6 1.64e-01 8.10e-02h 1
54r 0.0000000e+00 9.41e+00 5.06e+02 -4.2 5.20e+04 - 2.18e-02 8.75e-04h 1
55r 0.0000000e+00 8.46e+00 1.59e+03 -4.2 6.73e+03 - 2.61e-01 6.24e-02h 1
56r 0.0000000e+00 7.89e+00 1.55e+03 -4.2 9.99e+03 -8.1 4.90e-02 1.43e-02h 1
57r 0.0000000e+00 7.15e+00 1.47e+03 -4.2 2.18e+03 -7.7 4.69e-02 7.11e-02h 1
58r 0.0000000e+00 7.13e+00 1.47e+03 -4.2 6.27e+03 -8.1 3.08e-02 9.73e-04h 1
59r 0.0000000e+00 9.58e-01 5.60e+03 -4.2 3.41e+03 - 3.40e-01 5.07e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
60r 0.0000000e+00 9.23e-01 5.40e+03 -4.2 1.19e+03 - 1.00e+00 3.62e-02h 1
61r 0.0000000e+00 5.56e-01 1.77e+03 -4.2 7.65e+02 - 1.00e+00 6.99e-01h 1
62r 0.0000000e+00 6.45e-02 7.38e+01 -4.2 2.51e+02 - 1.00e+00 1.00e+00h 1
63r 0.0000000e+00 3.62e-03 3.12e-03 -4.2 3.02e+00 - 1.00e+00 1.00e+00h 1
64r 0.0000000e+00 3.62e-03 1.33e-08 -4.2 1.77e-03 - 1.00e+00 1.00e+00h 1
65r 0.0000000e+00 3.41e-04 6.48e-01 -6.4 7.54e+00 - 9.43e-01 9.64e-01f 1
66r 0.0000000e+00 1.31e+00 4.06e+02 -6.4 4.02e+03 -8.6 8.80e-01 2.16e-01f 1
67r 0.0000000e+00 1.12e+01 4.92e+02 -6.4 5.07e+03 -9.1 1.00e+00 4.71e-01h 1
68r 0.0000000e+00 1.10e+01 4.80e+02 -6.4 1.76e+04 - 1.81e-04 2.09e-02h 1
69r 0.0000000e+00 1.10e+01 4.79e+02 -6.4 1.88e+04 - 2.30e-04 1.15e-03h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
70r 0.0000000e+00 1.10e+01 4.79e+02 -6.4 1.45e+05 - 1.59e-07 6.05e-06f 1
71r 0.0000000e+00 1.10e+01 4.79e+02 -6.4 1.99e+04 - 4.84e-05 2.72e-05f 1
72r 0.0000000e+00 1.10e+01 4.79e+02 -6.4 8.49e+03 - 2.68e-01 9.54e-06f 1
73r 0.0000000e+00 1.00e+01 1.57e+03 -6.4 9.33e+03 - 8.71e-01 1.57e-01h 1
74r 0.0000000e+00 9.46e+00 8.10e+02 -6.4 5.23e+03 - 1.00e+00 7.50e-02h 1
75r 0.0000000e+00 9.04e+00 1.39e+03 -6.4 1.27e+04 - 4.86e-01 5.00e-02h 1
76r 0.0000000e+00 8.51e+00 6.82e+02 -6.4 2.05e+03 - 1.00e+00 5.84e-02h 1
77r 0.0000000e+00 1.54e+00 1.23e+02 -6.4 5.28e+02 - 1.00e+00 8.19e-01h 1
78r 0.0000000e+00 5.85e-05 4.22e+00 -6.4 1.50e+02 - 1.00e+00 1.00e+00h 1
79r 0.0000000e+00 5.85e-05 4.67e-04 -6.4 3.34e+00 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
80r 0.0000000e+00 5.85e-05 7.98e-10 -6.4 4.32e-03 - 1.00e+00 1.00e+00h 1
81r 0.0000000e+00 5.85e-05 1.74e+00 -9.0 2.08e-01 - 1.00e+00 9.59e-01f 1
82r 0.0000000e+00 1.96e-05 1.38e+02 -9.0 9.80e+02 - 1.00e+00 6.35e-01f 1
83r 0.0000000e+00 5.39e-09 1.67e+01 -9.0 3.58e+02 - 1.00e+00 8.79e-01f 1
Cannot recompute multipliers for feasibility problem. Error in eq_mult_calculator
Number of Iterations....: 83
(scaled) (unscaled)
Objective...............: 0.0000000000000000e+00 0.0000000000000000e+00
Dual infeasibility......: 0.0000000000000000e+00 0.0000000000000000e+00
Constraint violation....: 1.7678303265711293e-11 5.3885234052586384e-09
Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00
Overall NLP error.......: 1.7678303265711293e-11 5.3885234052586384e-09
Number of objective function evaluations = 86
Number of objective gradient evaluations = 5
Number of equality constraint evaluations = 96
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 85
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 83
Total CPU secs in IPOPT (w/o function evaluations) = 0.078
Total CPU secs in NLP function evaluations = 0.003
EXIT: Optimal Solution Found.
The costing block contains the custom model, and cost variables are in units of MUSD_2022 which match the default units of the custom class, self.base_currency = pyunits.USD_2022. Note that the model is not solved yet, so the variables are at their default values.
3.3 Build Process Costing#
Now, we can build the process costs. Because we defined the capital costing using a UnitModelCostingBlock, we see that the custom model is present in the registered cost for the flowsheet:
for block in m.fs.costing_2._registered_unit_costing:
print(block.name)
fs.custom_vessel.costing
This means that we don’t need to pass a special list including the unit, and the extra variables and constraints associated with the custom costing model will be added automatically. The REE Costing Framework will see that our custom model is already registered and will automatically check for capital_cost, fixed_operating_cost, and variable_operating_cost variables.
Note that we still need the product rates to calculate fixed operating costs.
Let’s build the process costs:
# Build process costs
m.fs.costing_2.build_process_costs(
Lang_factor=2.97,
fixed_OM=True,
pure_product_output_rates={
"treated_stream": m.fs.product[0] * 1.05
}, # includes injected water
mixed_product_output_rates={
"treated_stream": m.fs.product[0] * 0
}, # the product has 100% price realization, and is not considered a mixed product
sale_prices={"treated_stream": 1 * pyunits.USD_2021 / pyunits.gal},
variable_OM=True,
resources=[
"water",
"chemicals",
"waste",
],
rates=[m.fs.water, m.fs.chemicals, m.fs.waste],
prices={
"chemicals": 1.00 * pyunits.USD_2021 / pyunits.gallon,
"waste": 0.35 * pyunits.USD_2021 / pyunits.gallon,
},
watertap_blocks=[
m.fs.unit,
],
CE_index_year="2021",
)
QGESSCostingData.costing_initialization(m.fs.costing_2)
QGESSCostingData.initialize_fixed_OM_costs(m.fs.costing_2)
QGESSCostingData.initialize_variable_OM_costs(m.fs.costing_2)
# Solve and display results
solver.solve(m, tee=False)
# m.fs.costing_2.report() # uncomment to view results
assert m.fs.costing_2.total_BEC.value == pytest.approx(0.0038158, rel=1e-4)
assert m.fs.costing_2.total_installation_cost.value == pytest.approx(
0.0075172, rel=1e-4
)
assert m.fs.costing_2.total_plant_cost.value == pytest.approx(0.011333, rel=1e-4)
assert m.fs.costing_2.total_fixed_OM_cost.value == pytest.approx(5.1761, rel=1e-4)
assert m.fs.costing_2.total_sales_revenue.value == pytest.approx(0.085379, rel=1e-4)
assert m.fs.costing_2.total_variable_OM_cost[0].value == pytest.approx(3.9969, rel=1e-4)
assert m.fs.costing_2.plant_overhead_cost[0].value == pytest.approx(1.0352, rel=1e-4)
WARNING: model contains export suffix 'scaling_factor' that contains 27
component keys that are not exported as part of the NL file. Skipping.
WARNING: model contains export suffix 'scaling_factor' that contains 5 keys
that are not Var, Constraint, Objective, or the model. Skipping.
Now, let’s check some relevant costing model components to see how the custom model is integrated into the main costing equations:
for i in m.fs.costing_2.component_data_objects(Var, descend_into=True):
if "custom" in str(i):
print(i.name, "has a value of ", value(i))
print()
for i in m.fs.costing_2.component_data_objects(Constraint, descend_into=True):
if "custom" in str(i.expr):
print(i.expr)
print()
fs.costing_2.custom_fixed_costs has a value of 1.4700714402797727e-05
fs.costing_2.custom_variable_costs has a value of 8.495651334916658e-06
fs.costing_2.total_BEC == 1e-06*(MUSD_2021/USD_2021)*fs.custom_vessel.costing.capital_cost + 1.1739346708671863e-06*(MUSD_2021/USD_2018)*fs.unit.costing.capital_cost
fs.costing_2.custom_fixed_costs == 1.0000000000000002e-06*(MUSD_2021/USD_2021)*fs.custom_vessel.costing.fixed_operating_cost
fs.costing_2.total_fixed_OM_cost == fs.costing_2.annual_labor_cost + fs.costing_2.maintenance_and_material_cost + fs.costing_2.quality_assurance_and_control_cost + fs.costing_2.admin_and_support_labor_cost + fs.costing_2.sales_patenting_and_research_cost + fs.costing_2.property_taxes_and_insurance_cost + fs.costing_2.other_fixed_costs + fs.costing_2.watertap_fixed_costs + fs.costing_2.custom_fixed_costs
fs.costing_2.custom_variable_costs == 1.0000000000000002e-06*(MUSD_2021/USD_2021)*fs.custom_vessel.costing.variable_operating_cost
fs.costing_2.plant_overhead_cost[0.0] == 0.2*(fs.costing_2.total_fixed_OM_cost + (0*MUSD_2021)/a + (0*MUSD_2021)/a + (0*MUSD_2021)/a + fs.costing_2.custom_variable_costs)
fs.costing_2.total_variable_OM_cost[0.0] == fs.costing_2.variable_operating_costs[0.0,water] + fs.costing_2.variable_operating_costs[0.0,chemicals] + fs.costing_2.variable_operating_costs[0.0,waste] + fs.costing_2.other_variable_costs[0.0] + fs.costing_2.plant_overhead_cost[0.0] + (0*MUSD_2021)/a + (0*MUSD_2021)/a + (0*MUSD_2021)/a + fs.costing_2.custom_variable_costs
Although we didn’t change anything about the build_process_costs call to note the custom model, the costing framework extracted the cost variables from the custom vessel. The total_BEC expression is the capital cost of our custom vessel, since that is the only unit in the flowsheet, and we can see that it accounts for the unit conversion between the model units of MUSD_2022 and the flowsheet units of MUSD_2021. The operating cost equations also include the correct variables from our custom model, with the proper unit conversions.
4 Calculate Net Present Value and Taxes#
Now, let’s use the costing framework to calculate the net present value (NPV) of the process. The NPV represents the total value of the plant construction, operation, and production over the plant lifetime taking into account the changing worth of money over time. For example, distributed capital costs over multiple years or operating costs per year over the entire lifetime of a plant may grow with an escalation or inflation rate. The NPV provides a way to normalize costs over time to the current dollar year. Typical NPV components include capital cost, operating cost, revenue (this is a positive contribution, although we don’t have one in the current example), and cash flows from capital loan repayment.
To calculate the NPV as part of the build_process_costs method, certain configuration arguments must be set when the flowsheet costing block is instantiated. Users may provide their own values instead to bypass this requirement (see Section 4.1) or link to the costing framework (see Section 4.2). Additionally, the framework supports tax calculations from process revenue and costs (see Section 4.3).
4.1 Calculate NPV From Fixed Inputs#
The NPV calculation method may be utilized without calling the rest of the costing framework by setting fixed inputs for cost components when instantiating the costing block. This approach does not require any knowledge about the flowsheet, resources or product rates.
To begin, build the flowsheet costing block and specify the required cost components.
The total capital cost is 0.105 million USD, the total operating cost is the sum of the fixed and variable operating costs (9.173 million USD), and the total sales revenue is 0.085 million USD. The three arguments are required, and if one is set then all three must be set:
# Attach a flowsheet costing block with configuration arguments for NPV calculation
m.fs.costing_3 = QGESSCosting(
discount_percentage=10, # still a required argument, as before
plant_lifetime=20, # still a required argument, as before
# set values for capex, opex, and revenue
total_capital_cost=0.105, # capex
annual_operating_cost=9.173, # opex
annual_revenue=0.085, # revenue
cost_year="2021", # tells method to return results in MUSD_2021
# NPV arguments
has_capital_expenditure_period=True, # default is False, so need to set this to True
# use defaults for all other arguments, don't need to set them here
)
Note the total_capital_cost above is the entire capital cost including the total BEC, installation costs, and any other plant costs. This is not the same as the total_purchase_cost passed to the build_process_costs method, which comprises the total BEC only.
Users must set the argument cost_year to define the currency units. The arguments total_capital_cost, annual_operating_cost, and annual_revenue may be passed as scalars, but can also be passed as Expression, Var, or Param objects with or without units of measurement. If there are units of measurements, the variables will be converted internally to the units set by cost_year.
Let’s build the process costs. See below how we don’t need to define anything else about the plant, or even set any other arguments, to calculate the NPV:
m.fs.costing_3.build_process_costs(
fixed_OM=False, # the required components don't exist, so we skip the fixed_OM method
# by default, variable_OM is False
calculate_NPV=True,
CE_index_year="2021",
)
# Solve and display results (uncomment line to display)
solver.solve(m, tee=True)
# m.fs.costing_3.display() # uncomment to view results
assert m.fs.costing_3.npv.value == pytest.approx(-78.069, rel=1e-4)
WARNING: model contains export suffix 'scaling_factor' that contains 27
component keys that are not exported as part of the NL file. Skipping.
WARNING: model contains export suffix 'scaling_factor' that contains 5 keys
that are not Var, Constraint, Objective, or the model. Skipping.
Ipopt 3.13.2: nlp_scaling_method=gradient-based
tol=1e-06
max_iter=200
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt
This version of Ipopt was compiled from source code available at
https://github.com/IDAES/Ipopt as part of the Institute for the Design of
Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.
This version of Ipopt was compiled using HSL, a collection of Fortran codes
for large-scale scientific computation. All technical papers, sales and
publicity material resulting from use of the HSL codes within IPOPT must
contain the following acknowledgement:
HSL, a collection of Fortran codes for large-scale scientific
computation. See http://www.hsl.rl.ac.uk.
******************************************************************************
This is Ipopt version 3.13.2, running with linear solver ma27.
Number of nonzeros in equality constraint Jacobian...: 1273
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 529
Total number of variables............................: 437
variables with only lower bounds: 258
variables with lower and upper bounds: 146
variables with only upper bounds: 2
Total number of equality constraints.................: 437
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 0.0000000e+00 9.50e+03 1.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 0.0000000e+00 9.43e+03 3.03e+04 -1.0 4.01e+04 - 6.94e-03 7.56e-03h 1
2 0.0000000e+00 9.43e+03 3.05e+04 -1.0 3.92e+04 - 1.44e-02 7.65e-05h 1
3r 0.0000000e+00 9.43e+03 9.99e+02 2.0 0.00e+00 - 0.00e+00 3.83e-07R 2
4r 0.0000000e+00 8.91e+03 1.81e+04 2.0 2.30e+04 - 9.10e-03 6.61e-03f 1
5r 0.0000000e+00 8.35e+03 8.72e+03 2.0 6.40e+03 - 3.05e-02 8.06e-03f 1
6r 0.0000000e+00 7.83e+03 1.04e+04 2.0 4.79e+03 - 1.79e-01 1.02e-02f 1
7r 0.0000000e+00 2.51e+03 9.76e+03 2.0 3.14e+03 - 4.79e-01 2.39e-01f 1
8r 0.0000000e+00 7.39e+02 3.81e+03 2.0 1.84e+03 - 9.90e-01 9.07e-01f 1
9r 0.0000000e+00 4.48e+02 4.82e+03 1.3 5.24e+02 - 9.91e-01 6.56e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10r 0.0000000e+00 3.48e+02 5.56e+03 1.3 2.45e+03 - 1.00e+00 5.60e-01f 1
11r 0.0000000e+00 2.76e+02 1.41e+03 1.3 5.89e+02 - 1.00e+00 1.00e+00f 1
12r 0.0000000e+00 1.58e+02 6.35e+03 0.6 1.18e+03 - 9.40e-01 8.08e-01f 1
13r 0.0000000e+00 7.09e+01 3.52e+03 0.6 4.68e+02 - 6.94e-01 1.00e+00f 1
14r 0.0000000e+00 7.32e+01 5.22e+01 0.6 2.88e+02 - 1.00e+00 1.00e+00f 1
15r 0.0000000e+00 3.25e+01 2.86e+03 -0.8 1.85e+03 - 8.08e-01 7.12e-01f 1
16r 0.0000000e+00 2.57e+01 1.07e+04 -0.8 1.97e+02 - 9.75e-01 5.05e-01f 1
17r 0.0000000e+00 1.72e+01 4.03e+03 -0.8 4.67e+02 - 1.00e+00 7.65e-01f 1
18r 0.0000000e+00 1.34e+01 3.82e+03 -0.8 5.37e+02 - 1.00e+00 1.00e+00f 1
19r 0.0000000e+00 1.36e+01 1.03e+02 -0.8 1.01e+02 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20r 0.0000000e+00 1.36e+01 5.60e-02 -0.8 8.45e-02 - 1.00e+00 1.00e+00h 1
21r 0.0000000e+00 4.79e+00 3.75e+03 -3.3 1.24e+03 - 8.36e-01 6.73e-01f 1
22r 0.0000000e+00 2.84e-01 5.51e+03 -3.3 9.25e+02 - 5.07e-01 3.62e-01f 1
23r 0.0000000e+00 1.62e-01 7.81e+03 -3.3 8.16e+02 - 4.29e-02 2.10e-01f 1
24r 0.0000000e+00 1.35e-01 6.16e+03 -3.3 5.78e+02 - 1.94e-01 3.98e-01f 1
25r 0.0000000e+00 1.68e-01 2.70e+03 -3.3 3.20e+02 - 6.37e-01 4.78e-01f 1
26r 0.0000000e+00 1.11e-01 5.77e+02 -3.3 1.79e+02 - 5.77e-01 7.08e-01f 1
27r 0.0000000e+00 3.06e-02 2.52e+02 -3.3 4.37e+01 - 9.24e-01 7.99e-01f 1
28r 0.0000000e+00 4.62e-03 2.65e+01 -3.3 1.98e+01 - 1.00e+00 1.00e+00f 1
29r 0.0000000e+00 5.10e-03 5.80e-01 -3.3 3.67e+00 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
30r 0.0000000e+00 5.09e-03 5.95e-04 -3.3 2.16e-01 - 1.00e+00 1.00e+00h 1
31r 0.0000000e+00 5.51e-04 6.59e+00 -5.0 3.78e+00 - 9.23e-01 9.28e-01f 1
32r 0.0000000e+00 8.45e-03 3.12e+02 -5.0 1.09e+04 - 1.89e-02 3.58e-03f 1
33r 0.0000000e+00 1.15e-02 1.05e+03 -5.0 1.68e+03 - 6.50e-01 4.19e-02f 1
34r 0.0000000e+00 3.13e+00 1.01e+03 -5.0 2.17e+04 - 1.23e-02 6.04e-02f 1
35r 0.0000000e+00 3.07e+00 9.89e+02 -5.0 1.33e+01 -4.0 1.97e-01 1.99e-02h 1
36r 0.0000000e+00 2.85e+00 9.17e+02 -5.0 9.51e+00 -4.5 2.41e-01 7.32e-02h 1
37r 0.0000000e+00 2.25e+00 7.26e+02 -5.0 3.97e+01 -5.0 3.12e-02 2.10e-01h 1
38r 0.0000000e+00 1.88e+00 6.07e+02 -5.0 1.06e+02 -5.4 1.00e+00 1.63e-01h 1
39r 0.0000000e+00 1.86e+00 9.19e+02 -5.0 8.33e+03 - 3.41e-01 1.36e-02h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
40r 0.0000000e+00 1.84e+00 6.85e+02 -5.0 4.84e+03 - 1.37e-02 1.16e-01h 1
41r 0.0000000e+00 1.83e+00 1.84e+03 -5.0 3.85e+03 - 9.07e-01 9.95e-03h 1
42r 0.0000000e+00 1.82e+00 1.83e+03 -5.0 2.15e+04 - 2.85e-03 7.02e-03h 1
43r 0.0000000e+00 1.81e+00 1.78e+03 -5.0 5.94e+03 - 3.57e-02 3.18e-02h 2
44r 0.0000000e+00 1.82e+00 1.52e+03 -5.0 3.49e+03 - 3.92e-01 1.60e-01h 1
45r 0.0000000e+00 1.71e+00 1.42e+03 -5.0 1.09e+03 - 1.94e-01 6.15e-02h 1
46r 0.0000000e+00 1.66e+00 1.44e+03 -5.0 1.71e+03 - 8.54e-01 2.71e-02h 1
47r 0.0000000e+00 1.65e+00 1.31e+03 -5.0 3.48e+03 - 3.51e-01 8.69e-03h 1
48r 0.0000000e+00 1.47e+00 1.25e+03 -5.0 1.89e+03 - 6.46e-01 1.08e-01h 1
49r 0.0000000e+00 1.30e+00 4.98e+02 -5.0 1.56e+02 - 1.00e+00 1.12e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
50r 0.0000000e+00 1.34e-05 8.96e+00 -5.0 8.96e+01 - 1.00e+00 1.00e+00h 1
51r 0.0000000e+00 3.66e-06 4.88e-01 -5.0 4.85e+01 - 1.00e+00 1.00e+00h 1
52r 0.0000000e+00 1.23e-06 3.55e-03 -5.0 4.12e+00 - 1.00e+00 1.00e+00h 1
Cannot recompute multipliers for feasibility problem. Error in eq_mult_calculator
Number of Iterations....: 52
(scaled) (unscaled)
Objective...............: 0.0000000000000000e+00 0.0000000000000000e+00
Dual infeasibility......: 0.0000000000000000e+00 0.0000000000000000e+00
Constraint violation....: 2.9868630931691609e-08 1.2338951760337835e-06
Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00
Overall NLP error.......: 2.9868630931691609e-08 1.2338951760337835e-06
Number of objective function evaluations = 57
Number of objective gradient evaluations = 5
Number of equality constraint evaluations = 57
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 54
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 52
Total CPU secs in IPOPT (w/o function evaluations) = 0.044
Total CPU secs in NLP function evaluations = 0.004
EXIT: Optimal Solution Found.
As shown in the output above (uncomment display command to see output), only the NPV-related attributes are created.
Checking some cost results, we see that the costing block stored the cost inputs locally:
print(value(m.fs.costing_3.capex))
print(value(m.fs.costing_3.opex))
print(value(m.fs.costing_3.revenue))
0.105
9.173
0.085
The report method skips the attributes that don’t exist, so only the NPV is reported below:
m.fs.costing_3.report()
====================================================================================
costing_3
------------------------------------------------------------------------------------
Value
Plant Cost Units MUSD_2021
Net Present Value -78.069
====================================================================================
4.2 Calculate NPV From Costing Block#
The NPV calculations may also be coupled with the main costing framework, and automatically searches for relevant cost variables and expressions. To do this, we need to create our flowsheet costing block with several configuration arguments. Two arguments are required, the discount percentage defining cost projections over time and the plant lifetime in years. Many of the arguments have default values that are used if no other value is specified:
# Attach a flowsheet costing block with configuration arguments for NPV calculation
m.fs.costing_4 = QGESSCosting(
# arguments for NPV
discount_percentage=10, # percentage; Rate of return used to discount future cash flows back to their present value
plant_lifetime=20, # years; Length of operating period
has_capital_expenditure_period=True, # True/false flag whether a capital expenditure period occurs, default False
capital_expenditure_percentages=[
10,
60,
30,
], # A list of percentages that sum to 100 representing how capital costs are spread over a capital expenditure period
capital_escalation_percentage=3.6, # percentage; Rate at which capital costs escalate during the capital expenditure period, default 3.6%
capital_loan_interest_percentage=6, # percentage; Interest rate for capital equipment loan repayment, default 6%
capital_loan_repayment_period=10, # years; Length of loan repayment period, default 10 years
debt_percentage_of_capex=50, # percentage; Percentage of capex financed by debt, default 50%
operating_inflation_percentage=3, # percentage; Inflation rate for operating costs during the operating period, default 3%
revenue_inflation_percentage=3, # percentage; Inflation rate for revenue during the operating period, default 3%
)
If desired, users may set an argument debt_expression as an Expression to use instead of debt_percentage_of_capex * capex. This is an optional argument and will not be used here.
Now that the costing block has been created with the desired NPV-related arguments, we can set calculate_NPV in the build_process_costs method to add NPV calculations to the process costs. Because we did not specify any cost values in the costing block, the NPV method will look in the main costing block for required capital, operating, and revenue variables.
Now, let’s build the process costs, adding a new argument calculate_NPV set to True:
# Build process costs
m.fs.costing_4.build_process_costs(
Lang_factor=2.97,
fixed_OM=True,
# treated stream is not saleable, so pass dummy inputs
pure_product_output_rates={
"treated_stream": m.fs.product[0] * 1.05
}, # includes injected water
mixed_product_output_rates={
"treated_stream": m.fs.product[0] * 0
}, # the product has 100% price realization, and is not considered a mixed product
sale_prices={"treated_stream": 1 * pyunits.USD_2021 / pyunits.gal},
variable_OM=True,
resources=[
"water",
"chemicals",
"waste",
],
rates=[m.fs.water, m.fs.chemicals, m.fs.waste],
prices={
"chemicals": 1.00 * pyunits.USD_2021 / pyunits.gallon,
"waste": 0.35 * pyunits.USD_2021 / pyunits.gallon,
},
watertap_blocks=[
m.fs.unit,
],
calculate_NPV=True,
CE_index_year="2021",
)
QGESSCostingData.costing_initialization(m.fs.costing_4)
QGESSCostingData.initialize_fixed_OM_costs(m.fs.costing_4)
QGESSCostingData.initialize_variable_OM_costs(m.fs.costing_4)
# Solve and display results (uncomment line to display)
solver.solve(m, tee=True)
# m.fs.costing_4.display() # uncomment to view results
assert m.fs.costing_4.pv_capital_cost.value == pytest.approx(-0.00885464, rel=1e-4)
assert m.fs.costing_4.loan_debt.value == pytest.approx(0.00522987, rel=1e-4)
assert m.fs.costing_4.pv_loan_interest.value == pytest.approx(-0.000863716, rel=1e-4)
assert m.fs.costing_4.pv_operating_cost.value == pytest.approx(-78.7001, rel=1e-4)
assert m.fs.costing_4.pv_revenue.value == pytest.approx(0.7325064, rel=1e-4)
assert m.fs.costing_4.npv.value == pytest.approx(-77.9773, rel=1e-4)
WARNING: model contains export suffix 'scaling_factor' that contains 27
component keys that are not exported as part of the NL file. Skipping.
WARNING: model contains export suffix 'scaling_factor' that contains 5 keys
that are not Var, Constraint, Objective, or the model. Skipping.
Ipopt 3.13.2: nlp_scaling_method=gradient-based
tol=1e-06
max_iter=200
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt
This version of Ipopt was compiled from source code available at
https://github.com/IDAES/Ipopt as part of the Institute for the Design of
Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.
This version of Ipopt was compiled using HSL, a collection of Fortran codes
for large-scale scientific computation. All technical papers, sales and
publicity material resulting from use of the HSL codes within IPOPT must
contain the following acknowledgement:
HSL, a collection of Fortran codes for large-scale scientific
computation. See http://www.hsl.rl.ac.uk.
******************************************************************************
This is Ipopt version 3.13.2, running with linear solver ma27.
Number of nonzeros in equality constraint Jacobian...: 1341
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 529
Total number of variables............................: 464
variables with only lower bounds: 275
variables with lower and upper bounds: 148
variables with only upper bounds: 4
Total number of equality constraints.................: 464
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 0.0000000e+00 9.50e+03 1.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 0.0000000e+00 9.43e+03 3.03e+04 -1.0 4.01e+04 - 3.33e-03 7.56e-03h 1
2 0.0000000e+00 9.43e+03 3.05e+04 -1.0 3.92e+04 - 1.44e-02 7.65e-05h 1
3r 0.0000000e+00 9.43e+03 9.99e+02 2.3 0.00e+00 - 0.00e+00 3.83e-07R 2
4r 0.0000000e+00 9.35e+03 1.09e+04 2.3 3.74e+04 - 3.07e-03 2.14e-03f 1
5r 0.0000000e+00 9.21e+03 6.07e+03 2.3 2.91e+04 - 9.08e-03 4.03e-03f 1
6r 0.0000000e+00 9.03e+03 4.80e+03 2.3 5.88e+03 - 5.69e-02 5.37e-03f 1
7r 0.0000000e+00 8.10e+03 4.29e+03 2.3 3.59e+03 - 5.47e-02 3.42e-02f 1
8r 0.0000000e+00 6.05e+03 4.39e+03 2.3 2.88e+03 - 6.12e-02 1.03e-01f 1
9r 0.0000000e+00 5.09e+03 4.39e+03 2.3 2.20e+03 - 4.91e-02 8.81e-02f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10r 0.0000000e+00 4.27e+03 2.32e+04 2.3 1.84e+03 - 7.25e-01 1.14e-01f 1
11r 0.0000000e+00 1.97e+03 2.35e+04 2.3 1.11e+03 - 1.59e-01 1.00e+00f 1
12r 0.0000000e+00 2.61e+03 4.50e+02 2.3 1.27e+02 - 9.90e-01 1.00e+00f 1
13r 0.0000000e+00 7.46e+02 4.00e+03 0.9 6.94e+02 - 8.59e-01 7.48e-01f 1
14r 0.0000000e+00 6.12e+02 2.26e+04 0.9 5.31e+03 - 6.32e-01 1.86e-01f 1
15r 0.0000000e+00 3.05e+02 1.81e+04 0.9 3.42e+03 - 1.00e+00 6.21e-01f 1
16r 0.0000000e+00 1.29e+02 5.89e+03 0.9 7.66e+01 - 1.00e+00 9.31e-01f 1
17r 0.0000000e+00 1.30e+02 6.57e+02 0.9 5.16e+02 - 1.00e+00 1.00e+00f 1
18r 0.0000000e+00 5.37e+01 4.50e+03 0.2 1.67e+03 - 8.90e-01 8.45e-01f 1
19r 0.0000000e+00 4.10e+01 2.48e+02 0.2 8.73e+01 - 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20r 0.0000000e+00 4.20e+01 4.22e+01 0.2 2.50e+01 - 1.00e+00 1.00e+00f 1
21r 0.0000000e+00 2.31e+01 1.94e+03 -0.5 1.43e+03 - 1.00e+00 8.06e-01f 1
22r 0.0000000e+00 2.03e+01 5.96e+02 -0.5 5.25e+01 - 1.00e+00 1.00e+00f 1
23r 0.0000000e+00 2.02e+01 2.00e+02 -0.5 9.05e+01 - 1.00e+00 1.00e+00f 1
24r 0.0000000e+00 2.02e+01 5.13e+00 -0.5 5.62e+00 - 1.00e+00 1.00e+00h 1
25r 0.0000000e+00 1.06e+01 3.15e+03 -1.2 1.09e+03 - 1.00e+00 7.75e-01f 1
26r 0.0000000e+00 4.94e+00 1.79e+03 -1.2 3.02e+02 - 1.00e+00 1.00e+00f 1
27r 0.0000000e+00 4.60e+00 5.69e+01 -1.2 1.23e+02 - 1.00e+00 1.00e+00h 1
28r 0.0000000e+00 4.56e+00 5.08e-01 -1.2 2.24e-01 - 1.00e+00 1.00e+00h 1
29r 0.0000000e+00 2.24e-01 6.13e+03 -2.8 9.66e+02 - 9.26e-01 5.75e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
30r 0.0000000e+00 1.14e-01 6.28e+03 -2.8 5.24e+02 - 8.28e-01 3.65e-01f 1
31r 0.0000000e+00 1.23e-01 1.87e+03 -2.8 2.24e+02 - 1.00e+00 8.37e-01f 1
32r 0.0000000e+00 4.46e-02 2.14e+02 -2.8 5.84e+01 - 9.92e-01 1.00e+00f 1
33r 0.0000000e+00 3.74e-02 2.20e+01 -2.8 8.18e+00 - 1.00e+00 1.00e+00h 1
34r 0.0000000e+00 3.69e-02 4.35e-01 -2.8 1.07e+00 - 1.00e+00 1.00e+00h 1
35r 0.0000000e+00 3.69e-02 1.24e-04 -2.8 1.95e-02 - 1.00e+00 1.00e+00h 1
36r 0.0000000e+00 2.42e-02 4.39e+01 -6.4 3.10e+01 - 8.16e-01 8.55e-01f 1
37r 0.0000000e+00 8.56e-03 3.57e+02 -6.4 2.19e+01 -4.0 3.27e-01 1.88e-01f 1
38r 0.0000000e+00 7.45e-03 3.10e+02 -6.4 6.33e+01 -4.5 2.80e-01 1.31e-01f 1
39r 0.0000000e+00 6.89e-03 8.80e+02 -6.4 1.25e+02 -5.0 2.72e-01 7.70e-02f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
40r 0.0000000e+00 6.81e-03 1.39e+03 -6.4 3.01e+02 -5.4 1.85e-01 1.20e-02f 1
41r 0.0000000e+00 4.00e-02 6.39e+02 -6.4 6.27e+02 -5.9 1.41e-02 2.29e-01f 1
42r 0.0000000e+00 4.06e-02 6.34e+02 -6.4 1.54e+03 -6.4 1.40e-02 6.81e-03f 1
43r 0.0000000e+00 4.07e-02 6.33e+02 -6.4 2.19e+03 -6.9 1.37e-02 2.31e-03h 1
44r 0.0000000e+00 6.19e-02 1.02e+03 -6.4 1.38e+04 -7.3 8.27e-02 3.84e-03f 1
45r 0.0000000e+00 6.19e-02 1.04e+03 -6.4 2.52e+04 - 1.16e-01 8.43e-07h 1
46r 0.0000000e+00 6.90e-02 1.45e+03 -6.4 1.51e+04 - 1.51e-01 5.55e-03f 1
47r 0.0000000e+00 6.91e-02 1.62e+03 -6.4 6.82e+03 - 4.58e-01 4.86e-03h 1
48r 0.0000000e+00 7.44e-02 1.06e+03 -6.4 9.77e+02 - 7.68e-01 1.28e-01h 1
49r 0.0000000e+00 7.35e-02 1.35e+03 -6.4 3.90e+03 - 1.00e+00 1.27e-02h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
50r 0.0000000e+00 3.74e-02 4.42e+02 -6.4 1.67e+02 - 1.00e+00 4.92e-01h 1
51r 0.0000000e+00 1.59e-02 1.96e+02 -6.4 2.21e+02 - 1.00e+00 5.76e-01h 1
52r 0.0000000e+00 1.45e-03 7.10e+01 -6.4 1.01e+02 - 1.00e+00 9.09e-01f 1
53r 0.0000000e+00 5.61e-08 4.08e+01 -6.4 9.39e+00 - 9.15e-01 1.00e+00h 1
Cannot recompute multipliers for feasibility problem. Error in eq_mult_calculator
Number of Iterations....: 53
(scaled) (unscaled)
Objective...............: 0.0000000000000000e+00 0.0000000000000000e+00
Dual infeasibility......: 0.0000000000000000e+00 0.0000000000000000e+00
Constraint violation....: 3.1841977943258826e-09 5.6078025778560914e-08
Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00
Overall NLP error.......: 3.1841977943258826e-09 5.6078025778560914e-08
Number of objective function evaluations = 56
Number of objective gradient evaluations = 5
Number of equality constraint evaluations = 56
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 55
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 53
Total CPU secs in IPOPT (w/o function evaluations) = 0.046
Total CPU secs in NLP function evaluations = 0.000
EXIT: Optimal Solution Found.
As shown above (uncomment display command to see output), we’ve added several new variables to the model; present values are assumed to be negative unless otherwise noted. The attribute revenue is a Reference to total_sales_revenue and it is required because Pyomo bars assigning an object to two places in the same model.
The NPV method solves for the following variables:
pv_capital_cost- the present value of all capital costs over the plant lifetimepv_operating_cost- the present value of all operating costs over the plant lifetimepv_revenue- the present value of all revenue over the plant lifetime (positive value)loan_debt- principal loan taken out on capitalpv_loan_interest- the present value of all loan interest payments over the plant lifetimenpv- the net present value of the entire plant over its lifetime
and solves a new constraint corresponding to each new variable.
Let’s look at the costing report:
m.fs.costing_4.report()
====================================================================================
costing_4
------------------------------------------------------------------------------------
Value
Plant Cost Units MUSD_2021
Total Plant Cost 0.010460
Total Bare Erected Cost 0.0035218
Total Installation Cost 0.0069380
Total Other Plant Costs 1.0000e-12
Total Fixed Operating & Maintenance Cost 5.1760
Total Annual Operating Labor Cost 3.0730
Total Annual Technical Labor Cost 1.1801
Summation of Annual Labor Costs 4.2531
Total Maintenance and Material Cost 0.00020920
Total Quality Assurance and Control Cost 0.30730
Total Sales, Patenting and Research Cost 0.00042689
Summation of Sales, Admin and Insurance Cost 0.61513
Total Admin Support and Labor Cost 0.61460
Total Property Taxes and Insurance Cost 0.00010460
Total Other Fixed Costs 1.0000e-12
Total Variable Waste Cost 0.0000
Total Variable Chemicals Cost 0.0000
General Plant Overhead Cost 1.0352
Total Plant Overhead Cost, Including Maintenance & Quality Assurance 1.3427
Total Variable Operating & Maintenance Cost 3.9968
Total Land Cost 0.0000
Total Sales Revenue Cost 0.085377
Net Present Value -77.977
====================================================================================
We see that net present value is negative, which means that the plant is not economically viable and costs outweigh revenue over the plant lifetime.
Checking some cost results, we see that we see that that values very close to the prior case using fixed inputs, and that we obtain very similar NPV values:
print(value(m.fs.costing_4.capex))
print(value(m.fs.costing_4.opex))
print(value(m.fs.costing_4.revenue))
0.010459757918449586
9.172884943849951
0.08537719657954851
Case |
capex |
opex |
revenue |
NPV |
|---|---|---|---|---|
Fixed Input |
0.105 |
9.173 |
0.085 |
-78.069 |
Coupled With Costing |
0.010459757918449586 |
9.172884943849951 |
0.08537719657954851 |
-77.977 |
4.3 Tax and Incentives Calculations#
Suppose we know some information about the taxes and production incentives associated with constructing and operating this plant. We can define these properties, and the net present value method will automatically pick up these cash flows if present in the model. Users may include tax and incentives calculations by passing the following arguments to the build_process_costs method:
consider_taxes: True/False flag for calculating net tax owed. Defaults to False.income_tax_percentage: combined federal and state income tax percentage, usually between 26 - 40%. Here, it defaults to 26%.mineral_depletion_percentage: fixed tax deduction percentage for mineral depletion based on the type of mineral recovered, defaults to 14% of gross income excluding royalties as reported in the UKy report.production_incentive_percentage: tax deduction percentage for producing critical minerals, defaults to 10% of total production cost (excludes cost of feedstock).royalty_charge_percentage_of_revenue: Percentage of revenue charged as royalties; defaults to 6.5% as reported in the UKy report.
Income tax is a negative cash flow and is a percentage of the total net sales revenue (total gross sales revenue minus total annual production cost). Similarly, royalties are a negative cash flow and are a fixed percentage of the total gross sales revenue. Income tax is typically set by local laws applying to the construction site, whereas royalties are charged by the land owner. The default charge rates are 26% income tax, per NETL QGESS guidelines, and 6.5% royalties, per the University of Kentucky report.
Production incentives are a positive cash flow offsetting income tax based on the total annual production cost. Similarly, mineral depletion is a positive cash flow offsetting income tax based on total gross revenue excluding royalties. Production incentives are typically offered by governmental agencies to encourage production of certain material goods, whereas the mineral depletion charges allows plant operators to recover a portion of the expense associated with the reduced site value from extracting materials. The default rates are 10% production incentives, per the 45X mineral production credit, and 14% mineral depletion, per the University of Kentucky report.
Additionally, the net present value method will create and solve for a variable pv_taxes for the present value of taxes projected over the plant lifetime.
The University of Kentucky report and case study are discussed further in the UKy Flowsheet tutorial and with costing in the UKy Flowsheet with Costing tutorial.
As our current scenario involves brine nanofiltration and not mineral processing, the mineral depletion rate does not apply and is set to zero; however, the process still qualifies for a production incentive and this term is thus included.
The code below demonstrates including tax calculations in overnight cost and net present value calculations:
# Attach a new costing block to the flowsheet
m.fs.costing_5 = QGESSCosting(
# arguments for NPV
discount_percentage=10, # percent
plant_lifetime=20, # years
has_capital_expenditure_period=True,
capital_expenditure_percentages=[10, 60, 30],
capital_escalation_percentage=3.6,
capital_loan_interest_percentage=6,
capital_loan_repayment_period=10,
debt_percentage_of_capex=50,
operating_inflation_percentage=3,
revenue_inflation_percentage=3,
)
# Build process costs
m.fs.costing_5.build_process_costs(
Lang_factor=2.97,
fixed_OM=True,
# treated stream is not saleable, so pass dummy inputs
pure_product_output_rates={
"treated_stream": m.fs.product[0] * 1.05
}, # includes injected water
mixed_product_output_rates={
"treated_stream": m.fs.product[0] * 0
}, # the product has 100% price realization, and is not considered a mixed product
sale_prices={"treated_stream": 1 * pyunits.USD_2021 / pyunits.gal},
variable_OM=True,
resources=[
"water",
"chemicals",
"waste",
],
rates=[m.fs.water, m.fs.chemicals, m.fs.waste],
prices={
"chemicals": 1.00 * pyunits.USD_2021 / pyunits.gallon,
"waste": 0.35 * pyunits.USD_2021 / pyunits.gallon,
},
watertap_blocks=[
m.fs.unit,
],
consider_taxes=True,
income_tax_percentage=26,
mineral_depletion_percentage=0,
production_incentive_percentage=10,
royalty_charge_percentage_of_revenue=6.5,
calculate_NPV=True,
CE_index_year="2021",
)
QGESSCostingData.costing_initialization(m.fs.costing_5)
QGESSCostingData.initialize_fixed_OM_costs(m.fs.costing_5)
QGESSCostingData.initialize_variable_OM_costs(m.fs.costing_5)
solver.solve(m, tee=True)
# m.fs.costing_5.display() # uncomment to view results
assert value(m.fs.costing_5.royalty_charge) == pytest.approx(0.0055495, rel=1e-4)
assert value(m.fs.costing_5.mineral_depletion_charge) == pytest.approx(0, abs=1e-8)
assert value(m.fs.costing_5.production_incentive_charge) == pytest.approx(
0.91729, rel=1e-4
)
assert m.fs.costing_5.income_tax.value == pytest.approx(-2.3631, rel=1e-4)
assert m.fs.costing_5.additional_tax_owed.value == pytest.approx(0, abs=1e-8)
assert m.fs.costing_5.additional_tax_credit.value == pytest.approx(0, abs=1e-8)
assert m.fs.costing_5.net_tax_owed.value == pytest.approx(0, abs=1e-8)
WARNING: model contains export suffix 'scaling_factor' that contains 27
component keys that are not exported as part of the NL file. Skipping.
WARNING: model contains export suffix 'scaling_factor' that contains 5 keys
that are not Var, Constraint, Objective, or the model. Skipping.
Ipopt 3.13.2: nlp_scaling_method=gradient-based
tol=1e-06
max_iter=200
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt
This version of Ipopt was compiled from source code available at
https://github.com/IDAES/Ipopt as part of the Institute for the Design of
Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.
This version of Ipopt was compiled using HSL, a collection of Fortran codes
for large-scale scientific computation. All technical papers, sales and
publicity material resulting from use of the HSL codes within IPOPT must
contain the following acknowledgement:
HSL, a collection of Fortran codes for large-scale scientific
computation. See http://www.hsl.rl.ac.uk.
******************************************************************************
This is Ipopt version 3.13.2, running with linear solver ma27.
Number of nonzeros in equality constraint Jacobian...: 1423
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 544
Total number of variables............................: 494
variables with only lower bounds: 292
variables with lower and upper bounds: 150
variables with only upper bounds: 6
Total number of equality constraints.................: 494
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 0.0000000e+00 9.50e+03 1.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 0.0000000e+00 9.43e+03 3.03e+04 -1.0 4.01e+04 - 3.33e-03 7.56e-03h 1
2 0.0000000e+00 9.43e+03 3.05e+04 -1.0 3.92e+04 - 1.44e-02 7.65e-05h 1
3r 0.0000000e+00 9.43e+03 9.99e+02 2.4 0.00e+00 - 0.00e+00 3.83e-07R 2
4r 0.0000000e+00 9.38e+03 1.19e+04 2.4 3.76e+04 - 2.63e-03 1.70e-03f 1
5r 0.0000000e+00 9.24e+03 3.51e+03 2.4 3.27e+04 - 7.30e-03 4.86e-03f 1
6r 0.0000000e+00 9.14e+03 4.93e+03 2.4 6.08e+03 - 3.91e-02 3.39e-03f 1
7r 0.0000000e+00 8.47e+03 6.27e+03 2.4 3.84e+03 - 5.67e-02 2.67e-02f 1
8r 0.0000000e+00 6.77e+03 4.83e+03 2.4 2.88e+03 - 5.80e-02 8.94e-02f 1
9r 0.0000000e+00 5.74e+03 8.70e+03 2.4 2.18e+03 - 2.89e-01 8.85e-02f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10r 0.0000000e+00 4.57e+03 5.52e+03 2.4 1.52e+03 - 4.13e-02 2.12e-01f 1
11r 0.0000000e+00 2.71e+03 2.26e+04 2.4 1.13e+03 - 5.84e-02 5.81e-01f 1
12r 0.0000000e+00 2.81e+03 2.61e+04 2.4 3.97e+02 - 9.90e-01 3.25e-01f 1
13r 0.0000000e+00 8.08e+02 6.86e+03 1.7 6.67e+02 - 6.02e-01 8.97e-01f 1
14r 0.0000000e+00 5.27e+02 1.30e+04 1.0 2.09e+03 - 7.83e-01 3.99e-01f 1
15r 0.0000000e+00 3.59e+02 1.40e+04 1.0 3.76e+03 - 9.51e-01 6.03e-01f 1
16r 0.0000000e+00 1.35e+02 4.97e+03 1.0 1.49e+02 - 1.00e+00 9.76e-01f 1
17r 0.0000000e+00 1.46e+02 5.00e+02 1.0 4.55e+02 - 1.00e+00 1.00e+00f 1
18r 0.0000000e+00 6.28e+01 4.76e+03 0.3 1.57e+03 - 9.12e-01 8.48e-01f 1
19r 0.0000000e+00 4.40e+01 4.51e+02 0.3 1.33e+02 - 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20r 0.0000000e+00 4.51e+01 3.69e+01 0.3 5.88e+01 - 1.00e+00 1.00e+00f 1
21r 0.0000000e+00 2.42e+01 1.93e+03 -0.4 1.46e+03 - 1.00e+00 8.17e-01f 1
22r 0.0000000e+00 2.17e+01 4.89e+02 -0.4 7.07e+01 - 1.00e+00 1.00e+00f 1
23r 0.0000000e+00 2.17e+01 1.71e+02 -0.4 8.09e+01 - 1.00e+00 1.00e+00f 1
24r 0.0000000e+00 2.17e+01 4.38e+00 -0.4 5.83e+00 - 1.00e+00 1.00e+00h 1
25r 0.0000000e+00 1.22e+01 3.46e+03 -1.8 1.33e+03 - 8.42e-01 6.18e-01f 1
26r 0.0000000e+00 4.75e+00 8.85e+03 -1.8 6.63e+02 - 1.00e+00 6.33e-01f 1
27r 0.0000000e+00 2.74e-01 9.74e+03 -1.8 9.39e+02 - 1.00e+00 6.18e-01f 1
28r 0.0000000e+00 8.26e-02 2.87e+03 -1.8 3.16e+02 - 1.00e+00 7.65e-01f 1
29r 0.0000000e+00 5.75e-02 5.43e+01 -1.8 2.86e+01 - 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
30r 0.0000000e+00 5.69e-02 5.37e-01 -1.8 4.71e+00 - 1.00e+00 1.00e+00h 1
31r 0.0000000e+00 1.06e-01 2.02e+03 -2.7 1.48e+02 - 1.00e+00 5.88e-01f 1
32r 0.0000000e+00 8.50e-02 9.06e+02 -2.7 1.67e+02 - 7.28e-01 7.70e-01f 1
33r 0.0000000e+00 4.59e-02 4.43e+01 -2.7 2.95e+01 - 1.00e+00 1.00e+00f 1
34r 0.0000000e+00 4.50e-02 4.32e-01 -2.7 1.06e+00 - 1.00e+00 1.00e+00h 1
35r 0.0000000e+00 4.50e-02 3.71e-04 -2.7 3.68e-02 - 1.00e+00 1.00e+00h 1
36r 0.0000000e+00 3.39e-02 4.86e+01 -6.2 3.36e+01 - 8.12e-01 8.26e-01f 1
37r 0.0000000e+00 1.10e-02 3.67e+02 -6.2 2.27e+01 -4.0 3.32e-01 1.86e-01f 1
38r 0.0000000e+00 8.95e-03 2.93e+02 -6.2 6.69e+01 -4.5 2.54e-01 2.01e-01f 1
39r 0.0000000e+00 8.55e-03 9.22e+02 -6.2 1.26e+02 -5.0 2.90e-01 4.41e-02f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
40r 0.0000000e+00 8.14e-03 1.30e+03 -6.2 2.75e+02 -5.4 1.99e-01 5.89e-02f 1
41r 0.0000000e+00 4.05e-02 6.65e+02 -6.2 5.23e+02 -5.9 4.45e-02 2.49e-01f 1
42r 0.0000000e+00 4.13e-02 6.55e+02 -6.2 1.07e+03 -6.4 2.31e-02 1.24e-02f 1
43r 0.0000000e+00 4.13e-02 6.51e+02 -6.2 1.07e+03 -6.9 3.80e-02 4.62e-03h 1
44r 0.0000000e+00 6.05e-02 1.00e+03 -6.2 1.84e+03 -7.3 5.29e-01 2.87e-02f 1
45r 0.0000000e+00 6.21e-02 1.27e+03 -6.2 1.23e+04 - 1.27e-01 3.71e-03h 1
46r 0.0000000e+00 6.19e-02 1.36e+03 -6.2 5.95e+03 - 6.15e-01 8.05e-03h 1
47r 0.0000000e+00 6.58e-02 1.17e+03 -6.2 3.21e+03 - 8.72e-01 8.69e-02h 1
48r 0.0000000e+00 5.69e-02 1.20e+03 -6.2 3.07e+03 - 1.00e+00 1.48e-01h 1
49r 0.0000000e+00 3.60e-02 5.23e+02 -6.2 4.61e+02 - 1.00e+00 3.69e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
50r 0.0000000e+00 8.76e-03 1.37e+02 -6.2 3.10e+02 - 1.00e+00 7.57e-01h 1
51r 0.0000000e+00 8.33e-04 8.94e+01 -6.2 9.28e+01 - 1.00e+00 9.05e-01f 1
52r 0.0000000e+00 1.24e-07 2.68e+01 -6.2 9.16e+00 - 9.43e-01 1.00e+00f 1
Cannot recompute multipliers for feasibility problem. Error in eq_mult_calculator
Number of Iterations....: 52
(scaled) (unscaled)
Objective...............: 0.0000000000000000e+00 0.0000000000000000e+00
Dual infeasibility......: 0.0000000000000000e+00 0.0000000000000000e+00
Constraint violation....: 1.0220808643879309e-09 1.2443839558784475e-07
Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00
Overall NLP error.......: 1.0220808643879309e-09 1.2443839558784475e-07
Number of objective function evaluations = 55
Number of objective gradient evaluations = 5
Number of equality constraint evaluations = 55
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 54
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 52
Total CPU secs in IPOPT (w/o function evaluations) = 0.056
Total CPU secs in NLP function evaluations = 0.001
EXIT: Optimal Solution Found.
m.fs.costing_5.report()
====================================================================================
costing_5
------------------------------------------------------------------------------------
Value
Plant Cost Units MUSD_2021
Total Plant Cost 0.010460
Total Bare Erected Cost 0.0035218
Total Installation Cost 0.0069380
Total Other Plant Costs 1.0000e-12
Total Fixed Operating & Maintenance Cost 5.1760
Total Annual Operating Labor Cost 3.0730
Total Annual Technical Labor Cost 1.1801
Summation of Annual Labor Costs 4.2531
Total Maintenance and Material Cost 0.00020920
Total Quality Assurance and Control Cost 0.30730
Total Sales, Patenting and Research Cost 0.00042689
Summation of Sales, Admin and Insurance Cost 0.61513
Total Admin Support and Labor Cost 0.61460
Total Property Taxes and Insurance Cost 0.00010460
Total Other Fixed Costs 1.0000e-12
Total Variable Waste Cost 0.0000
Total Variable Chemicals Cost 0.0000
General Plant Overhead Cost 1.0352
Total Plant Overhead Cost, Including Maintenance & Quality Assurance 1.3427
Total Variable Operating & Maintenance Cost 3.9968
Total Land Cost 0.0000
Total Sales Revenue Cost 0.085377
Net Present Value -77.977
Royalty Charge 0.0055495
Mineral Depletion Charge 0.0000
Production Incentive Charge 0.91741
Income Tax -2.3631
Additional Tax Owed 1.0000e-12
Additional Tax Credit 1.0000e-12
Net Tax Owed 7.6439e-10
====================================================================================
In the report above, the tax components are now included. Processes in which the production costs are large and the sales revenue is small could have a large production credit and a small combined income tax and royalties cash flow, yielding a negative tax. In the case of a calculated negative tax, the tax contribution is bounded by an internal model variable min_net_tax_owed that defaults to 0. Users may override this enforce a minimum tax owed, or provide an allowance for a negative tax. The minimum tax value is changed below and the resulting net tax owed and NPV are presented:
print("Net Taxed Owed: ", value(m.fs.costing_5.net_tax_owed))
print("NPV: ", value(m.fs.costing_5.npv))
print("\nEnforce a minimum tax of $3 million per year")
m.fs.costing_5.min_net_tax_owed.fix(3)
solver.solve(m, tee=False)
print("Net Taxed Owed: ", value(m.fs.costing_5.net_tax_owed))
print("NPV: ", value(m.fs.costing_5.npv))
print(
"\nArtificially increase the tax credit and allow a negative tax up to -$1 million per year"
)
m.fs.costing_5.additional_tax_credit.fix(5)
m.fs.costing_5.min_net_tax_owed.fix(-1)
solver.solve(m, tee=False)
print("Net Taxed Owed: ", value(m.fs.costing_5.net_tax_owed))
print("NPV: ", value(m.fs.costing_5.npv))
Net Taxed Owed: 7.643892704437963e-10
NPV: -77.97734966988045
Enforce a minimum tax of $3 million per year
WARNING: model contains export suffix 'scaling_factor' that contains 27
component keys that are not exported as part of the NL file. Skipping.
WARNING: model contains export suffix 'scaling_factor' that contains 5 keys
that are not Var, Constraint, Objective, or the model. Skipping.
Net Taxed Owed: 3.000000000396878
NPV: -97.16645263036031
Artificially increase the tax credit and allow a negative tax up to -$1 million per year
WARNING: model contains export suffix 'scaling_factor' that contains 27
component keys that are not exported as part of the NL file. Skipping.
WARNING: model contains export suffix 'scaling_factor' that contains 5 keys
that are not Var, Constraint, Objective, or the model. Skipping.
Net Taxed Owed: -0.9999999996557521
NPV: -71.58102814994122
5 Economy of Numbers#
The costing framework supports calculation of process cost savings using economy of numbers, which estimates future profitability of novel equipment and cost savings due to experiential learning. Let’s suppose that our brine nanofiltration process becomes more efficient over time as the technology develops; we want to see how much the cost reduces due to knowledge gain over time.
5.1 Mathematical Definition#
The cost of manufacturing a piece of equipment tends to decline as the production quantity increases (economy of scales), but also as new generations of that equipment are developed leading to increases in technical knowledge (economy of numbers). New equipment are often denoted as First-of-a-Kind (FOAK) equipment while later generations with improved knowledge are denoted as Nth-of-a-Kind (NOAK) equipment.
The following equations are taken from Faber G, Ruttinger A, Strunge T, Langhorst T, Zimmermann A, van der Hulst M, Bensebaa F, Moni S and Tao L (2022) Adapting Technology Learning Curves for Prospective Techno-Economic and Life Cycle Assessments of Emerging Carbon Capture and Utilization Pathways. Front. Clim. 4:820261. doi: 10.3389/fclim.2022.820261:
\(Y = A * x^b\)
\(b = - \frac{\log_{10}(1-R)}{\log_{10}(2)}\)
where \(Y\) is the NOAK cost, \(A\) is the FOAK cost, \(x\) is the cumulative number of units (or generations of equipment), \(b\) is the learning rate exponent, and \(R\) is the learning rate. The learning rate is a matter of engineering knowledge for various technologies and plant configurations, and range between 0.01 to 0.1 depending on the level of maturity (experimental, growing, proven, etc.). More information on typical learning rates may be found at Rubin, E. S., Mantripragada, H., and Zhai, H., “An Assessment of the NETL Cost Estimation Methodology”. Department of Engineering and Public Policy, Carnegie Mellon University, Pittsburgh, PA (2016). p. 31, Fig. 6-4.
5.2 Usage Example#
Our first-generation brine nanofiltration resulted in a total (installed) plant cost of $0.10459 million. The item has undergone development and is a 5th-generation design. The learning rate for the technology is 0.05.
The following code may be used to predict the NOAK cost:
# Calculate the economy of numbers benefit
QGESSCostingData.economy_of_numbers(
m.fs.costing_5,
cum_num_units=5,
cost_FOAK=m.fs.costing_5.total_plant_cost,
CE_index_year="2021",
learning_rate=0.05,
)
solver = get_solver()
results = solver.solve(m, tee=True)
WARNING: model contains export suffix 'scaling_factor' that contains 27
component keys that are not exported as part of the NL file. Skipping.
WARNING: model contains export suffix 'scaling_factor' that contains 5 keys
that are not Var, Constraint, Objective, or the model. Skipping.
Ipopt 3.13.2: nlp_scaling_method=gradient-based
tol=1e-06
max_iter=200
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt
This version of Ipopt was compiled from source code available at
https://github.com/IDAES/Ipopt as part of the Institute for the Design of
Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.
This version of Ipopt was compiled using HSL, a collection of Fortran codes
for large-scale scientific computation. All technical papers, sales and
publicity material resulting from use of the HSL codes within IPOPT must
contain the following acknowledgement:
HSL, a collection of Fortran codes for large-scale scientific
computation. See http://www.hsl.rl.ac.uk.
******************************************************************************
This is Ipopt version 3.13.2, running with linear solver ma27.
Number of nonzeros in equality constraint Jacobian...: 1425
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 544
Total number of variables............................: 495
variables with only lower bounds: 293
variables with lower and upper bounds: 150
variables with only upper bounds: 6
Total number of equality constraints.................: 495
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 0.0000000e+00 1.00e+05 1.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 0.0000000e+00 9.92e+04 3.03e+04 -1.0 1.00e+05 - 6.94e-03 7.56e-03h 1
2 0.0000000e+00 9.92e+04 3.05e+04 -1.0 9.92e+04 - 1.44e-02 7.65e-05h 1
3r 0.0000000e+00 9.92e+04 9.99e+02 5.0 0.00e+00 - 0.00e+00 3.83e-07R 2
4r 0.0000000e+00 9.92e+04 4.71e+04 5.0 3.33e+05 - 2.03e-03 2.15e-06f 1
5r 0.0000000e+00 8.30e+04 4.88e+04 3.6 3.13e+05 - 9.24e-01 5.20e-02f 1
6r 0.0000000e+00 1.74e+04 4.58e+04 1.5 1.40e+06 - 1.21e-01 5.88e-02f 1
7r 0.0000000e+00 1.71e+04 4.51e+04 1.5 4.75e+04 - 6.13e-03 1.73e-02f 1
8r 0.0000000e+00 1.69e+04 4.43e+04 1.5 6.96e+02 - 1.93e-01 1.43e-02f 1
9r 0.0000000e+00 1.16e+04 3.01e+04 1.5 2.66e+02 - 1.74e-02 3.28e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10r 0.0000000e+00 9.74e+03 2.49e+04 1.5 1.80e+02 - 1.93e-01 1.72e-01f 1
11r 0.0000000e+00 8.28e+03 2.13e+04 1.5 2.22e+02 - 4.79e-02 1.63e-01f 1
12r 0.0000000e+00 7.07e+03 1.84e+04 1.5 1.76e+02 - 1.31e-01 1.60e-01f 1
13r 0.0000000e+00 6.03e+03 1.63e+04 1.5 1.06e+02 - 1.06e-01 1.66e-01f 1
14r 0.0000000e+00 4.24e+03 1.95e+04 1.5 8.89e+01 - 1.33e-01 3.41e-01f 1
15r 0.0000000e+00 2.60e+03 2.49e+04 1.5 6.98e+01 - 4.39e-01 4.76e-01f 1
16r 0.0000000e+00 1.27e+03 3.39e+04 1.5 5.70e+01 - 2.81e-01 7.27e-01f 1
17r 0.0000000e+00 8.84e+02 1.13e+04 1.5 1.03e+02 - 5.95e-01 7.24e-01f 1
18r 0.0000000e+00 5.92e+02 6.06e+03 0.8 4.89e+02 - 5.44e-01 3.87e-01f 1
19r 0.0000000e+00 4.30e+02 7.03e+04 0.8 1.21e+03 - 8.16e-02 7.48e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20r 0.0000000e+00 3.75e+02 5.20e+04 0.8 1.93e+02 - 4.58e-01 2.91e-01f 1
21r 0.0000000e+00 1.30e+02 2.45e+04 0.8 2.33e+02 - 4.44e-01 7.82e-01f 1
22r 0.0000000e+00 1.03e+02 3.06e+04 0.8 3.83e+02 - 6.06e-01 9.00e-01f 1
23r 0.0000000e+00 1.04e+02 1.04e+04 0.8 5.04e+02 - 6.50e-01 1.00e+00f 1
24r 0.0000000e+00 1.05e+02 4.85e+03 0.8 2.06e+03 - 1.00e+00 1.00e+00f 1
25r 0.0000000e+00 1.05e+02 7.77e+02 0.8 8.38e+02 - 1.00e+00 1.00e+00f 1
26r 0.0000000e+00 1.05e+02 1.06e+00 0.8 1.27e+00 - 1.00e+00 1.00e+00h 1
27r 0.0000000e+00 4.09e+01 3.65e+03 -2.0 1.94e+03 - 7.55e-01 6.85e-01f 1
28r 0.0000000e+00 3.05e+01 8.74e+03 -2.0 3.28e+02 - 7.80e-01 4.92e-01f 1
29r 0.0000000e+00 2.30e+01 8.56e+03 -2.0 2.80e+02 - 7.14e-01 4.50e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
30r 0.0000000e+00 1.41e+01 6.14e+03 -2.0 8.10e+02 - 7.74e-01 4.85e-01f 1
31r 0.0000000e+00 4.97e+00 1.08e+04 -2.0 1.70e+03 - 2.99e-01 6.15e-01f 1
32r 0.0000000e+00 1.51e+00 6.99e+03 -2.0 1.52e+03 - 8.01e-01 4.93e-01f 1
33r 0.0000000e+00 1.77e-01 9.64e+03 -2.0 9.82e+02 - 9.06e-01 6.04e-01f 1
34r 0.0000000e+00 8.51e-02 4.89e+03 -2.0 2.08e+02 - 1.00e+00 6.98e-01f 1
35r 0.0000000e+00 3.78e-02 2.88e+02 -2.0 4.71e+01 - 1.00e+00 1.00e+00f 1
36r 0.0000000e+00 3.66e-02 3.62e+00 -2.0 6.68e+00 - 1.00e+00 1.00e+00h 1
37r 0.0000000e+00 3.66e-02 1.17e-03 -2.0 9.45e-02 - 1.00e+00 1.00e+00h 1
38r 0.0000000e+00 9.40e-02 1.41e+03 -4.5 1.14e+02 - 9.87e-01 6.35e-01f 1
39r 0.0000000e+00 2.59e-02 1.35e+03 -4.5 1.64e+03 - 9.33e-02 1.44e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
40r 0.0000000e+00 2.05e-02 1.63e+03 -4.5 7.33e+02 - 9.51e-01 8.26e-02f 1
41r 0.0000000e+00 8.78e-02 1.07e+03 -4.5 7.68e+02 - 2.89e-01 4.35e-01f 1
42r 0.0000000e+00 5.90e-02 7.28e+02 -4.5 4.11e+02 - 1.00e+00 4.17e-01f 1
43r 0.0000000e+00 2.41e-02 4.67e+02 -4.5 2.43e+02 - 1.00e+00 1.00e+00f 1
44r 0.0000000e+00 4.06e-05 3.16e+00 -4.5 1.63e+01 - 1.00e+00 1.00e+00h 1
45r 0.0000000e+00 3.19e-05 9.64e-02 -4.5 2.77e+00 - 1.00e+00 1.00e+00h 1
46r 0.0000000e+00 3.20e-05 1.26e-04 -4.5 1.02e-01 - 1.00e+00 1.00e+00h 1
47r 0.0000000e+00 3.29e-05 1.05e+00 -6.7 6.91e+00 - 9.20e-01 9.33e-01f 1
48r 0.0000000e+00 3.24e-05 4.82e+01 -6.7 4.29e+00 -4.0 1.11e-01 1.05e-01f 1
49r 0.0000000e+00 2.23e-04 1.06e+03 -6.7 5.58e+01 -4.5 2.06e-02 2.35e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
50r 0.0000000e+00 2.03e-04 1.03e+03 -6.7 1.56e+02 -5.0 1.17e-03 2.92e-02f 1
51r 0.0000000e+00 2.09e-04 1.02e+03 -6.7 4.18e+02 -5.4 1.86e-02 5.84e-03f 1
52r 0.0000000e+00 2.09e-04 1.02e+03 -6.7 2.56e+02 -5.9 1.92e-02 1.51e-03h 1
53r 0.0000000e+00 2.14e-04 1.94e+03 -6.7 6.52e+02 -6.4 4.53e-01 1.26e-02f 1
54r 0.0000000e+00 2.12e-04 1.91e+03 -6.7 3.68e+02 -6.9 1.51e-02 1.07e-02f 1
55r 0.0000000e+00 6.05e-04 1.91e+03 -6.7 1.41e+05 - 4.37e-04 5.08e-04f 1
56r 0.0000000e+00 6.05e-04 1.91e+03 -6.7 6.50e+02 -7.3 4.35e-03 7.25e-04f 1
57r 0.0000000e+00 7.74e-04 1.89e+03 -6.7 3.39e+03 -7.8 4.52e-02 2.04e-03f 1
58r 0.0000000e+00 7.69e-04 1.56e+03 -6.7 7.93e+03 - 9.68e-01 1.61e-02h 1
59r 0.0000000e+00 8.81e-04 1.45e+03 -6.7 7.27e+03 - 1.00e+00 3.82e-02h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
60r 0.0000000e+00 7.74e-04 1.65e+03 -6.7 3.83e+02 - 1.00e+00 1.23e-01h 1
61r 0.0000000e+00 1.88e-07 1.41e+00 -6.7 4.94e+00 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 61
(scaled) (unscaled)
Objective...............: 0.0000000000000000e+00 0.0000000000000000e+00
Dual infeasibility......: 0.0000000000000000e+00 0.0000000000000000e+00
Constraint violation....: 1.2509190128184855e-08 1.8762932718630054e-07
Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00
Overall NLP error.......: 1.2509190128184855e-08 1.8762932718630054e-07
Number of objective function evaluations = 64
Number of objective gradient evaluations = 5
Number of equality constraint evaluations = 64
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 63
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 61
Total CPU secs in IPOPT (w/o function evaluations) = 0.055
Total CPU secs in NLP function evaluations = 0.003
EXIT: Optimal Solution Found.
# display results (uncomment line to display)
m.fs.costing_5.display()
assert m.fs.costing_5.cost_NOAK.value == pytest.approx(0.00928533, rel=1e-4)
Block fs.costing_5
Variables:
total_BEC : Total TPC
Size=1, Index=None, Units=MUSD_2021
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : 0 : 0.003521804012602152 : None : False : False : Reals
total_installation_cost : Total installation cost
Size=1, Index=None, Units=MUSD_2021
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : 0 : 0.006937953904826436 : None : False : False : Reals
total_plant_cost : Total plant cost
Size=1, Index=None, Units=MUSD_2021
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : 0 : 0.010459757918428758 : None : False : False : Reals
other_plant_costs : Additional plant costs
Size=1, Index=None, Units=MUSD_2021
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : 0 : 1e-12 : None : True : True : Reals
annual_operating_labor_cost : Annual operating labor cost
Size=1, Index=None, Units=MUSD_2021/a
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : 0 : 3.072988800000001 : None : False : False : Reals
annual_technical_labor_cost : Annual technical labor cost
Size=1, Index=None, Units=MUSD_2021/a
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : 0 : 1.1800656 : None : False : False : Reals
annual_labor_cost : Annual labor cost
Size=1, Index=None, Units=MUSD_2021/a
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : 0 : 4.253054400000002 : None : False : False : Reals
maintenance_and_material_cost : Maintenance and material cost
Size=1, Index=None, Units=MUSD_2021/a
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : 0 : 0.00020919515836891343 : None : False : False : Reals
quality_assurance_and_control_cost : Quality assurance and control cost
Size=1, Index=None, Units=MUSD_2021/a
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : 0 : 0.3072988800000001 : None : False : False : Reals
sales_patenting_and_research_cost : Sales, patenting and research cost
Size=1, Index=None, Units=MUSD_2021/a
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : 0 : 0.0004268860592769847 : None : False : False : Reals
admin_and_support_labor_cost : Admin and support labor cost
Size=1, Index=None, Units=MUSD_2021/a
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : 0 : 0.6145977600000002 : None : False : False : Reals
property_taxes_and_insurance_cost : Property taxes and insurance cost
Size=1, Index=None, Units=MUSD_2021/a
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : 0 : 0.00010459757918496274 : None : False : False : Reals
total_fixed_OM_cost : Total fixed O&M costs
Size=1, Index=None, Units=MUSD_2021/a
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : 0 : 5.176043899207736 : None : False : False : Reals
total_sales_revenue : Total sales revenue
Size=1, Index=None, Units=MUSD_2021/a
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : 0 : 0.08537721185536334 : None : False : False : Reals
other_fixed_costs : Other fixed costs
Size=1, Index=None, Units=MUSD_2021/a
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : 0 : 1e-12 : None : True : True : Reals
watertap_fixed_costs : WaterTAP fixed costs
Size=1, Index=None, Units=MUSD_2021/a
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : 0 : 0.0003521804012603573 : None : False : False : Reals
custom_fixed_costs : Custom fixed costs
Size=1, Index=None, Units=MUSD_2021/a
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : 0 : 8.642507838324042e-12 : None : False : False : Reals
variable_operating_costs : Variable operating costs
Size=3, Index=fs._time*{water, chemicals, waste}, Units=MUSD_2021/a
Key : Lower : Value : Upper : Fixed : Stale : Domain
(0.0, 'chemicals') : None : 0.16128000000000003 : None : False : False : Reals
(0.0, 'waste') : None : 2.7692252202601355 : None : False : False : Reals
(0.0, 'water') : None : 0.031127040000000005 : None : False : False : Reals
other_variable_costs : A variable to include non-standard O&M costs
Size=1, Index=fs._time, Units=MUSD_2021/a
Key : Lower : Value : Upper : Fixed : Stale : Domain
0.0 : 0 : 1e-12 : None : True : True : Reals
custom_variable_costs : Custom variable costs
Size=1, Index=None, Units=MUSD_2021/a
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : 0 : 8.642507588803356e-12 : None : False : False : Reals
total_variable_OM_cost : Total variable operating and maintenance costs
Size=1, Index=fs._time, Units=MUSD_2021/a
Key : Lower : Value : Upper : Fixed : Stale : Domain
0.0 : None : 3.996841040113054 : None : False : False : Reals
plant_overhead_cost : Plant overhead costs
Size=1, Index=fs._time, Units=MUSD_2021/a
Key : Lower : Value : Upper : Fixed : Stale : Domain
0.0 : None : 1.0352087798432756 : None : False : False : Reals
net_tax_owed : Net tax owed in millions USD
Size=1, Index=None, Units=MUSD_2021/a
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : None : -0.9999999996563536 : None : False : False : Reals
income_tax : Income tax in millions USD
Size=1, Index=None, Units=MUSD_2021/a
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : None : -2.3630637464106155 : None : False : False : Reals
additional_tax_credit : Additional tax credit
Size=1, Index=None, Units=MUSD_2021/a
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : None : 5 : None : True : True : Reals
additional_tax_owed : Additional tax owed
Size=1, Index=None, Units=MUSD_2021/a
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : None : 1e-12 : None : True : True : Reals
min_net_tax_owed : Minimum net tax owed in millions USD
Size=1, Index=None, Units=MUSD_2021/a
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : None : -1 : None : True : True : Reals
pv_capital_cost : Present value of total lifetime capital costs; negative cash flow
Size=1, Index=None, Units=MUSD_2021
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : None : -0.008854635374064565 : 0 : False : False : Reals
loan_debt : total debt from loans in $MM
Size=1, Index=None, Units=MUSD_2021
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : 0 : 0.005229878959214308 : 10000.0 : False : False : Reals
pv_loan_interest : present value of total lifetime loan interest in $MM; normally a negative cash flow, but can be positive depending on the discount and interest rates
Size=1, Index=None, Units=MUSD_2021
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : -10000.0 : -0.0008637156251232432 : 10000.0 : False : False : Reals
pv_operating_cost : Present value of total lifetime operating costs; negative cash flow
Size=1, Index=None, Units=MUSD_2021
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : None : -78.70013737323399 : 0 : False : False : Reals
pv_revenue : Present value of total lifetime sales revenue; positive cash flow
Size=1, Index=None, Units=MUSD_2021
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : 0 : 0.732506550121223 : None : False : False : Reals
pv_taxes : Present value of total lifetime tax owed; typically negative cash flow, positive in the case of an effective negative tax rate
Size=1, Index=None, Units=MUSD_2021
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : None : 6.396366428875211 : None : False : False : Reals
npv : Present value of plant over entire capital and operation lifetime
Size=1, Index=None, Units=MUSD_2021
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : None : -71.58098274523675 : None : False : False : Reals
cost_NOAK : Cost of the Nth-of-A-Kind of the unit
Size=1, Index=None, Units=MUSD_2021
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : 0 : 0.00928533221110632 : None : False : False : Reals
Objectives:
None
Constraints:
total_BEC_eq : Size=1
Key : Lower : Body : Upper
None : 0.0 : 5.928417479150738e-16 : 0.0
total_installation_cost_eq : Size=1
Key : Lower : Body : Upper
None : 0.0 : 1.960237527853792e-16 : 0.0
total_plant_cost_eq : Size=1
Key : Lower : Body : Upper
None : 0.0 : 1.700029006457271e-16 : 0.0
annual_operating_labor_cost_eq : Size=1
Key : Lower : Body : Upper
None : 3.072988800000001 : 3.072988800000001 : 3.072988800000001
annual_technical_labor_cost_eq : Size=1
Key : Lower : Body : Upper
None : 1.1800656 : 1.1800656 : 1.1800656
annual_labor_cost_eq : Size=1
Key : Lower : Body : Upper
None : 0.0 : 8.881784197001252e-16 : 0.0
maintenance_and_material_cost_eq : Size=1
Key : Lower : Body : Upper
None : 0.0 : 3.382710778154774e-16 : 0.0
quality_assurance_and_control_cost_eq : Size=1
Key : Lower : Body : Upper
None : 0.0 : 0.0 : 0.0
sales_patenting_and_research_cost_eq : Size=1
Key : Lower : Body : Upper
None : 0.0 : 1.679971266266289e-16 : 0.0
admin_and_support_labor_cost_eq : Size=1
Key : Lower : Body : Upper
None : 0.0 : 0.0 : 0.0
taxes_and_insurance_cost_eq : Size=1
Key : Lower : Body : Upper
None : 0.0 : 6.751597978610357e-16 : 0.0
sum_watertap_fixed_costs : Size=1
Key : Lower : Body : Upper
None : 0.0 : 2.0144476364780672e-16 : 0.0
sum_custom_fixed_costs : Size=1
Key : Lower : Body : Upper
None : 1e-12 : 8.642507838324042e-12 : 1e-12
total_fixed_OM_cost_eq : Size=1
Key : Lower : Body : Upper
None : 0.0 : 0.0 : 0.0
total_sales_revenue_eq : Size=1
Key : Lower : Body : Upper
None : 0.0 : 4.163336342344337e-17 : 0.0
variable_cost_eq : Size=3
Key : Lower : Body : Upper
(0.0, 'chemicals') : 0.0 : 0.0 : 0.0
(0.0, 'waste') : 0.0 : 0.0 : 0.0
(0.0, 'water') : 0.0 : 0.0 : 0.0
sum_custom_variable_costs : Size=1
Key : Lower : Body : Upper
None : 1e-12 : 8.642507588803356e-12 : 1e-12
plant_overhead_cost_eq : Size=1
Key : Lower : Body : Upper
0.0 : 0.0 : -2.220446049250313e-16 : 0.0
total_variable_cost_eq : Size=1
Key : Lower : Body : Upper
0.0 : 0.0 : 4.440892098500626e-16 : 0.0
income_tax_eq : Size=1
Key : Lower : Body : Upper
None : 0.0 : 4.440892098500626e-16 : 0.0
net_tax_owed_eq : Size=1
Key : Lower : Body : Upper
None : 0.0 : 0.0 : 0.0
pv_capital_cost_constraint : Size=1
Key : Lower : Body : Upper
None : 0.0 : -8.673617379884035e-18 : 0.0
loan_debt_constraint : Size=1
Key : Lower : Body : Upper
None : 0.0 : 1.3877787807814457e-17 : 0.0
pv_loan_interest_constraint : Size=1
Key : Lower : Body : Upper
None : 0.0 : 1.0842021724855044e-19 : 0.0
pv_operating_cost_constraint : Size=1
Key : Lower : Body : Upper
None : 0.0 : 0.0 : 0.0
pv_revenue_constraint : Size=1
Key : Lower : Body : Upper
None : 0.0 : 0.0 : 0.0
pv_taxes_constraint : Size=1
Key : Lower : Body : Upper
None : 0.0 : 0.0 : 0.0
npv_constraint : Size=1
Key : Lower : Body : Upper
None : 0.0 : -1.4210854715202004e-14 : 0.0
cost_NOAK_eq : Size=1
Key : Lower : Body : Upper
None : 0.0 : 1.700029006457271e-16 : 0.0
As shown in the output above (uncomment display command to see output), the 5th-generation equipment cost is $ 0.09285 million USD, which corresponds to a cost reduction of 11.2% relative to the first-generation cost of $ 0.10459 million USD.
6 Costing Bounding#
Finally, the costing framework supports calculating costing bounds based on the feedstock capacity and grade of a process. This method applies to various types of mineral processing, not brine nanofiltration; as such, a toy example will be used here rather than the brine nanofiltration process
Like the fixed input NPV method in section 4 and the economy of numbers method in Section 5, the costing bounding method is independent of the main costing framework. The bound equations derive from Fritz, A.G., Tarka, T.J. & Mauter, M.S. Assessing the economic viability of unconventional rare earth element feedstocks. Nat Sustain 6, 1103–1112 (2023). https://doi.org/10.1038/s41893-023-01145-1. The study establishes a database of capital and operating expenses for REE recovery and market prices of saleable mineral products, and develops surrogates models to predict the maximum cost threshold of specific process types to determine the economic viability. Thresholds are calculated using a power law:
\( C_{threshold} (2022 USD / kg REE) = \alpha * x^\beta\)
where \(C_{threshold}\) is the cost threshold, \(x\) is a known quantity about the process, and \(\alpha\) and \(\beta\) are empirical coefficients. A dictionary of coefficients were regressed from process data for a range of process types with \(x\) as the total feedstock REE:
\( x = F_{grade}(unitless) * F_{capacity}(tonnes)\)
The values [\(\alpha\), \(\beta\)] are available in the public paper and the costing framework code, and are presented below:
“Total Capital”: [81, -0.46]
“Total Operating”: [27, -0.087]
“Beneficiation”: [2.7, -0.15]
“Beneficiation, Chemical Extraction, Enrichment and Separation”: [22, -0.059]
“Chemical Extraction”: [40, -0.46]
“Chemical Extraction, Enrichment and Separation”: [15, -0.19]
“Enrichment and Separation”: [6.7, -0.16]
“Mining”: [25, -0.32]
6.1 Quantifying Threshold Bounds#
To quantify uncertainty in the threshold estimation, a 95% confidence interval was applied to the values to produce upper and lower bound predictions for the cost threshold per the following equation:
\(C_{threshold}^{\pm} (2022 USD / kg REE) = ({\alpha} \pm {\alpha}_{0.95}) * x^{\beta \pm {\beta}_{0.95}}\)
where \(C_{threshold}^{-}\) corresponds to the 5% confidence bound and \(C_{threshold}^{+}\) corresponds to the 95% confidence bound.
The following values, available in the costing framework code, were produced for [\(\alpha\), \({\alpha}_{0.95}\), \(\beta\), \({\beta}_{0.95}\)]:
“Total Capital”: [81, 1.4, -0.46, 0.063]
“Total Operating”: [27, 0.87, -0.087, 0.038]
“Beneficiation”: [2.7, 1.3, -0.15, 0.062]
“Beneficiation, Chemical Extraction, Enrichment and Separation”: [22, 1.28, -0.059, 0.046]
“Chemical Extraction, Enrichment and Separation”: [15, 15, -0.19, 0.28]
“Enrichment and Separation”: [6.7, 2.8, -0.16, 0.11]
“Mining”: [25, 2.5, -0.32, 0.095]
6.2 Usage Example#
The method itself is straightforward to implement if the required inputs are provided. Suppose we have a process with a feed input rate of 500 U.S. ton of REE per hour with an REE grade of 356.64 ppm. We can call the costing bounding method, which will print the calculated bound values after solving:
# Create a Concrete Model as the top level object
# m = ConcreteModel()
# Add a flowsheet object to the model
# m.fs = FlowsheetBlock(dynamic=False)
# Attach a flowsheet costing block with configuration arguments for NPV calculation
m.fs.costing_6 = QGESSCosting()
# Define inputs
m.fs.feed_input = Var(initialize=500, units=pyunits.ton / pyunits.hr)
m.fs.feed_grade = Var(initialize=356.64, units=pyunits.ppm)
# Define annual operating hours
hours_per_shift = 8
shifts_per_day = 3
operating_days_per_year = 336
m.fs.annual_operating_hours = Param(
initialize=hours_per_shift * shifts_per_day * operating_days_per_year,
mutable=True,
units=pyunits.hours / pyunits.a, # pyunits.a is "annum", equivalent to pyunits.year
)
# Define plant lifetime
m.fs.plant_lifetime = Param(initialize=20, mutable=True, units=pyunits.a)
# Call costing bound method
CE_index_year = "2021"
QGESSCostingData.calculate_REE_costing_bounds(
b=m.fs.costing_6, # block to attach costing bound variables and constraints to
capacity=m.fs.feed_input
* m.fs.annual_operating_hours
* m.fs.plant_lifetime, # lifetime capacity
grade=m.fs.feed_grade, # feed grade
CE_index_year=CE_index_year, # project dollar year
)
2025-08-19 08:02:04 [INFO] idaes.prommis.uky.costing.ree_plant_capcost: New variable 'capacity' created as attribute of fs.costing_6
2025-08-19 08:02:04 [INFO] idaes.prommis.uky.costing.ree_plant_capcost: New variable 'grade' created as attribute of fs.costing_6
2025-08-19 08:02:04 [INFO] idaes.prommis.uky.costing.ree_plant_capcost: New variable 'costing_lower_bound' created as attribute of fs.costing_6
2025-08-19 08:02:04 [INFO] idaes.prommis.uky.costing.ree_plant_capcost: New variable 'costing_upper_bound' created as attribute of fs.costing_6
2025-08-19 08:02:04 [INFO] idaes.prommis.uky.costing.ree_plant_capcost: New constraint 'costing_lower_bounding_eq' created as attribute of fs.costing_6
2025-08-19 08:02:04 [INFO] idaes.prommis.uky.costing.ree_plant_capcost: New constraint 'costing_upper_bounding_eq' created as attribute of fs.costing_6
2025-08-19 08:02:04 [INFO] idaes.prommis.uky.costing.ree_plant_capcost:
Printing calculated costing bounds for processes:
Total Capital : [ 0.3384067822001138 , 1.2616218937057968 ] USD_2021 /kg
Total Operating : [ 6.359499980920717 , 14.691722903144084 ] USD_2021 /kg
Beneficiation : [ 0.14066244356767998 , 1.4182512509318754 ] USD_2021 /kg
Beneficiation, Chemical Extraction, Enrichment and Separation : [ 6.180204384847299 , 17.69749844860086 ] USD_2021 /kg
Chemical Extraction : [ 0.07208272971633867 , 1.4372356232012462 ] USD_2021 /kg
Chemical Extraction, Enrichment and Separation : [ 0.0 , 65.00506864477617 ] USD_2021 /kg
Enrichment and Separation : [ 0.21724965888728387 , 4.957273599035815 ] USD_2021 /kg
Mining : [ 0.28687186039287327 , 2.420854142560117 ] USD_2021 /kg
The results above were calculated directly from the inputs, and in this case are the correct values. This is not the case for flowsheets in general, as we did not solve the flowsheet yet.
Solving and printing the bound values, we confirm the same results:
# Solve and display results
solver.solve(m, tee=True)
processes = [
"Total Capital",
"Total Operating",
"Beneficiation",
"Beneficiation, Chemical Extraction, Enrichment and Separation",
"Chemical Extraction",
"Chemical Extraction, Enrichment and Separation",
"Enrichment and Separation",
"Mining",
]
print("\n\nPrinting calculated costing bounds for processes:")
for p in processes:
print(
p,
": [",
value(m.fs.costing_6.costing_lower_bound[p]),
", ",
value(m.fs.costing_6.costing_upper_bound[p]),
"]",
getattr(pyunits, "USD_" + CE_index_year),
"/kg",
)
# check a few entries
assert m.fs.costing_6.costing_lower_bound["Total Capital"].value == pytest.approx(
0.338407, rel=1e-4
)
assert m.fs.costing_6.costing_upper_bound["Total Capital"].value == pytest.approx(
1.26162, rel=1e-4
)
assert m.fs.costing_6.costing_lower_bound["Total Operating"].value == pytest.approx(
6.35950, rel=1e-4
)
assert m.fs.costing_6.costing_upper_bound["Total Operating"].value == pytest.approx(
14.6917, rel=1e-4
)
WARNING: model contains export suffix 'scaling_factor' that contains 27
component keys that are not exported as part of the NL file. Skipping.
WARNING: model contains export suffix 'scaling_factor' that contains 5 keys
that are not Var, Constraint, Objective, or the model. Skipping.
Ipopt 3.13.2: nlp_scaling_method=gradient-based
tol=1e-06
max_iter=200
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt
This version of Ipopt was compiled from source code available at
https://github.com/IDAES/Ipopt as part of the Institute for the Design of
Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.
This version of Ipopt was compiled using HSL, a collection of Fortran codes
for large-scale scientific computation. All technical papers, sales and
publicity material resulting from use of the HSL codes within IPOPT must
contain the following acknowledgement:
HSL, a collection of Fortran codes for large-scale scientific
computation. See http://www.hsl.rl.ac.uk.
******************************************************************************
This is Ipopt version 3.13.2, running with linear solver ma27.
Number of nonzeros in equality constraint Jacobian...: 1441
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 544
Total number of variables............................: 511
variables with only lower bounds: 309
variables with lower and upper bounds: 150
variables with only upper bounds: 6
Total number of equality constraints.................: 511
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 0.0000000e+00 9.50e+03 1.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 0.0000000e+00 9.43e+03 3.04e+04 -1.0 4.01e+04 - 6.94e-03 7.56e-03h 1
2 0.0000000e+00 9.43e+03 3.05e+04 -1.0 3.92e+04 - 1.44e-02 7.65e-05h 1
3r 0.0000000e+00 9.43e+03 9.99e+02 0.3 0.00e+00 - 0.00e+00 3.83e-07R 2
4r 0.0000000e+00 7.92e+02 9.97e+02 0.3 5.81e+04 - 1.51e-03 2.03e-03f 1
5r 0.0000000e+00 3.74e+02 9.97e+02 0.3 3.32e+03 - 5.39e-04 1.98e-03f 1
6r 0.0000000e+00 1.72e+02 3.53e+03 0.3 1.36e+04 - 2.19e-02 2.36e-03f 1
7r 0.0000000e+00 9.46e+01 6.19e+03 0.3 2.03e+04 - 2.70e-02 1.36e-02f 1
8r 0.0000000e+00 6.45e+01 1.46e+04 0.3 1.28e+04 - 1.01e-01 2.60e-02f 1
9r 0.0000000e+00 4.89e+01 2.18e+04 0.3 8.28e+03 - 1.58e-01 8.16e-02f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10r 0.0000000e+00 4.33e+01 2.33e+04 0.3 3.62e+03 - 3.56e-01 1.55e-01f 1
11r 0.0000000e+00 4.28e+01 4.61e+04 0.3 1.72e+03 - 8.23e-01 1.71e-01f 1
12r 0.0000000e+00 4.28e+01 3.44e+04 0.3 1.06e+03 - 8.92e-01 3.69e-01f 1
13r 0.0000000e+00 4.67e+01 1.17e+04 0.3 9.70e+02 - 9.90e-01 8.84e-01f 1
14r 0.0000000e+00 4.80e+01 4.29e+03 0.3 7.55e+01 - 4.86e-01 8.96e-01f 1
15r 0.0000000e+00 4.67e+01 4.38e+02 0.3 5.12e+01 - 1.00e+00 1.00e+00f 1
16r 0.0000000e+00 4.69e+01 1.87e+00 0.3 6.93e+00 - 1.00e+00 1.00e+00f 1
17r 0.0000000e+00 2.51e+01 2.80e+03 -1.8 1.82e+03 - 8.65e-01 6.53e-01f 1
18r 0.0000000e+00 1.89e+01 6.16e+03 -1.8 1.13e+02 - 8.65e-01 5.19e-01f 1
19r 0.0000000e+00 1.23e+01 1.18e+04 -1.8 8.14e+02 - 9.90e-01 4.42e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20r 0.0000000e+00 2.57e+00 8.87e+03 -1.8 1.51e+03 - 1.00e+00 6.96e-01f 1
21r 0.0000000e+00 2.37e-01 9.23e+03 -1.8 9.60e+02 - 9.57e-01 6.27e-01f 1
22r 0.0000000e+00 2.15e-01 2.70e+03 -1.8 1.33e+02 - 1.00e+00 8.18e-01f 1
23r 0.0000000e+00 6.26e-02 4.32e+01 -1.8 2.39e+01 - 1.00e+00 1.00e+00f 1
24r 0.0000000e+00 6.28e-02 3.53e-01 -1.8 9.51e-01 - 1.00e+00 1.00e+00h 1
25r 0.0000000e+00 1.06e-01 2.27e+03 -2.7 1.63e+02 - 1.00e+00 5.60e-01f 1
26r 0.0000000e+00 8.09e-02 1.01e+03 -2.7 1.64e+02 - 7.68e-01 7.84e-01f 1
27r 0.0000000e+00 4.84e-02 3.37e+01 -2.7 2.61e+01 - 1.00e+00 1.00e+00f 1
28r 0.0000000e+00 4.73e-02 2.49e-01 -2.7 5.52e-01 - 1.00e+00 1.00e+00h 1
29r 0.0000000e+00 4.73e-02 3.75e-05 -2.7 3.88e-02 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
30r 0.0000000e+00 3.75e-02 6.48e+01 -6.1 3.49e+01 - 8.14e-01 8.08e-01f 1
31r 0.0000000e+00 1.20e-02 3.55e+02 -6.1 2.29e+01 -4.0 3.38e-01 1.95e-01f 1
32r 0.0000000e+00 9.41e-03 2.72e+02 -6.1 6.77e+01 -4.5 2.47e-01 2.32e-01f 1
33r 0.0000000e+00 9.03e-03 9.16e+02 -6.1 1.26e+02 -5.0 3.03e-01 4.03e-02f 1
34r 0.0000000e+00 8.50e-03 1.28e+03 -6.1 2.61e+02 -5.4 2.14e-01 7.60e-02f 1
35r 0.0000000e+00 4.08e-02 6.79e+02 -6.1 4.76e+02 -5.9 4.32e-02 2.63e-01f 1
36r 0.0000000e+00 4.16e-02 6.67e+02 -6.1 8.89e+02 -6.4 2.92e-02 1.62e-02f 1
37r 0.0000000e+00 4.15e-02 6.60e+02 -6.1 7.18e+02 -6.9 1.12e-01 6.07e-03h 1
38r 0.0000000e+00 5.03e-02 6.54e+02 -6.1 1.50e+04 - 8.15e-02 5.66e-03f 1
39r 0.0000000e+00 5.06e-02 1.03e+03 -6.1 1.09e+04 - 2.09e-01 1.75e-03f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
40r 0.0000000e+00 5.05e-02 1.31e+03 -6.1 5.32e+03 - 6.33e-01 2.37e-03f 1
41r 0.0000000e+00 6.45e-02 1.24e+03 -6.1 4.28e+03 - 9.44e-01 1.14e-01h 1
42r 0.0000000e+00 6.15e-02 1.20e+03 -6.1 2.19e+03 - 1.00e+00 4.96e-02h 1
43r 0.0000000e+00 3.84e-02 6.43e+02 -6.1 1.26e+03 - 1.00e+00 3.92e-01h 1
44r 0.0000000e+00 7.05e-03 1.03e+02 -6.1 1.89e+02 - 1.00e+00 8.22e-01h 1
45r 0.0000000e+00 6.79e-04 9.57e+01 -6.1 9.88e+01 - 1.00e+00 9.05e-01f 1
46r 0.0000000e+00 9.49e-08 1.75e+01 -6.1 1.08e+01 - 9.62e-01 1.00e+00f 1
Number of Iterations....: 46
(scaled) (unscaled)
Objective...............: 0.0000000000000000e+00 0.0000000000000000e+00
Dual infeasibility......: 0.0000000000000000e+00 0.0000000000000000e+00
Constraint violation....: 7.7696142852801131e-08 9.4895066027333996e-08
Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00
Overall NLP error.......: 7.7696142852801131e-08 9.4895066027333996e-08
Number of objective function evaluations = 49
Number of objective gradient evaluations = 5
Number of equality constraint evaluations = 49
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 48
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 46
Total CPU secs in IPOPT (w/o function evaluations) = 0.044
Total CPU secs in NLP function evaluations = 0.001
EXIT: Optimal Solution Found.
Printing calculated costing bounds for processes:
Total Capital : [ 0.33840678220011383 , 1.2616218937057968 ] USD_2021 /kg
Total Operating : [ 6.359499980920717 , 14.691722903144084 ] USD_2021 /kg
Beneficiation : [ 0.14066244356768 , 1.4182512509318754 ] USD_2021 /kg
Beneficiation, Chemical Extraction, Enrichment and Separation : [ 6.180204384847299 , 17.69749844860086 ] USD_2021 /kg
Chemical Extraction : [ 0.07208272971633865 , 1.4372356232012462 ] USD_2021 /kg
Chemical Extraction, Enrichment and Separation : [ 3.353678230577125e-10 , 65.00506864477617 ] USD_2021 /kg
Enrichment and Separation : [ 0.21724965888728387 , 4.957273599035815 ] USD_2021 /kg
Mining : [ 0.2868718603928733 , 2.420854142560117 ] USD_2021 /kg
Summary#
The advanced method demonstrated here expand upon the capabilities of the REE Costing Framework. Syntax for the example WaterTAP and custom models may be applied to any imported WaterTAP or custom model, respectively, and it is recommended that additional unit operation specificaitions are placed prior to calling the build_process_costs. The NPV calculations must be set up when instantiating the flowsheet costing block itself, as shown in this notebook, while the costing bounding equations are entirely independent of the capital and operating cost equations and may be added at any point while building the flowsheet.