Hda Flowsheet With Costing¶
HDA Flowsheet Costing¶
Note¶
This tutorial will demonstrate adding capital and operating costs to the two HDA examples, the basic HDA with Flash and a comparison with the HDA with Distillation.
Learning outcomes¶
- Import external pre-built steady-state flowsheets using the IDAES unit model library
- Define and add costing blocks using the IDAES Process Costing Framework
- Fomulate and solve a process economics optimization problem
- Defining an objective function
- Setting variable bounds
- Adding additional constraints
Problem Statement¶
Hydrodealkylation is a chemical reaction that often involves reacting an aromatic hydrocarbon in the presence of hydrogen gas to form a simpler aromatic hydrocarbon devoid of functional groups. In this example, toluene will be reacted with hydrogen gas at high temperatures to form benzene via the following reaction:
C6H5CH3 + H2 → C6H6 + CH4
This reaction is often accompanied by an equilibrium side reaction which forms diphenyl, which we will neglect for this example.
This example is based on the 1967 AIChE Student Contest problem as present by Douglas, J.M., Chemical Design of Chemical Processes, 1988, McGraw-Hill.
Users may refer to the prior examples linked at the top of this notebook for detailed process descriptions of the two HDA configurations. As before, the properties required for this module are defined in
hda_ideal_VLE.py
idaes.models.properties.activity_coeff_models.BTX_activity_coeff_VLE
hda_reaction.py
Additionally, we will be importing externally-defined flowsheets for the two HDA configurations from
hda_flowsheets_for_costing_notebook.py
Import and run HDA Flowsheets¶
First, we will generate solved flowsheets for each HDA model. The external scripts build and set inputs for the flowsheets, initialize unit models and streams, and solve the flowsheets before returning the model objects. Note that the HDA flowsheets contain all unit models and stream connections, and no costing equations.
The flowsheet utilizes the Wegstein method to iteratively solve circular dependencies such as recycle streams, and is intended to approach a feasible solution. As such, the calls below will fail to converge after 3 iterations and pass to IPOPT to obtain an optimal solution as expected:
# Source file for prebuilt flowsheets
from hda_flowsheets_for_costing_notebook import hda_with_flash
# Build hda model with second flash unit and return model object
m = hda_with_flash(tee=True)
Building flowsheet... Setting inputs... Initializing flowsheet... Limiting Wegstein tear to 3 iterations to obtain initial solution, if not converged IPOPT will pick up and continue. fs.s03 fs.H101 fs.R101 fs.F101 fs.S101 fs.C101 fs.M101 2023-03-04 01:43:25 [INFO] idaes.init.fs.H101.control_volume: Initialization Complete 2023-03-04 01:43:25 [INFO] idaes.init.fs.H101: Initialization Complete: optimal - Optimal Solution Found 2023-03-04 01:43:25 [INFO] idaes.init.fs.R101.control_volume: Initialization Complete 2023-03-04 01:43:25 [INFO] idaes.init.fs.R101: Initialization Complete: optimal - Optimal Solution Found 2023-03-04 01:43:25 [INFO] idaes.init.fs.F101.control_volume: Initialization Complete 2023-03-04 01:43:25 [INFO] idaes.init.fs.F101: Initialization Complete: optimal - Optimal Solution Found 2023-03-04 01:43:25 [INFO] idaes.init.fs.S101.purge_state: Initialization Complete 2023-03-04 01:43:25 [INFO] idaes.init.fs.S101.recycle_state: Initialization Complete 2023-03-04 01:43:25 [INFO] idaes.init.fs.S101: Initialization Step 2 Complete: optimal - Optimal Solution Found 2023-03-04 01:43:25 [INFO] idaes.init.fs.F102.control_volume: Initialization Complete 2023-03-04 01:43:25 [INFO] idaes.init.fs.F102: Initialization Complete: optimal - Optimal Solution Found 2023-03-04 01:43:25 [INFO] idaes.init.fs.C101.control_volume: Initialization Complete 2023-03-04 01:43:25 [INFO] idaes.init.fs.C101: Initialization Complete: optimal - Optimal Solution Found 2023-03-04 01:43:25 [INFO] idaes.init.fs.M101.mixed_state: Initialization Complete 2023-03-04 01:43:25 [INFO] idaes.init.fs.M101: Initialization Complete: optimal - Optimal Solution Found 2023-03-04 01:43:25 [INFO] idaes.init.fs.H101.control_volume: Initialization Complete 2023-03-04 01:43:25 [INFO] idaes.init.fs.H101: Initialization Complete: optimal - Optimal Solution Found 2023-03-04 01:43:25 [INFO] idaes.init.fs.R101.control_volume: Initialization Complete 2023-03-04 01:43:25 [INFO] idaes.init.fs.R101: Initialization Complete: optimal - Optimal Solution Found 2023-03-04 01:43:25 [INFO] idaes.init.fs.F101.control_volume: Initialization Complete 2023-03-04 01:43:26 [INFO] idaes.init.fs.F101: Initialization Complete: optimal - Optimal Solution Found 2023-03-04 01:43:26 [INFO] idaes.init.fs.S101.purge_state: Initialization Complete 2023-03-04 01:43:26 [INFO] idaes.init.fs.S101.recycle_state: Initialization Complete 2023-03-04 01:43:26 [INFO] idaes.init.fs.S101: Initialization Step 2 Complete: optimal - Optimal Solution Found 2023-03-04 01:43:26 [INFO] idaes.init.fs.C101.control_volume: Initialization Complete 2023-03-04 01:43:26 [INFO] idaes.init.fs.C101: Initialization Complete: optimal - Optimal Solution Found 2023-03-04 01:43:26 [INFO] idaes.init.fs.M101.mixed_state: Initialization Complete 2023-03-04 01:43:26 [INFO] idaes.init.fs.M101: Initialization Complete: optimal - Optimal Solution Found 2023-03-04 01:43:26 [INFO] idaes.init.fs.H101.control_volume: Initialization Complete 2023-03-04 01:43:26 [INFO] idaes.init.fs.H101: Initialization Complete: optimal - Optimal Solution Found 2023-03-04 01:43:26 [INFO] idaes.init.fs.R101.control_volume: Initialization Complete 2023-03-04 01:43:26 [INFO] idaes.init.fs.R101: Initialization Complete: optimal - Optimal Solution Found 2023-03-04 01:43:26 [INFO] idaes.init.fs.F101.control_volume: Initialization Complete 2023-03-04 01:43:26 [INFO] idaes.init.fs.F101: Initialization Complete: optimal - Optimal Solution Found 2023-03-04 01:43:26 [INFO] idaes.init.fs.S101.purge_state: Initialization Complete 2023-03-04 01:43:26 [INFO] idaes.init.fs.S101.recycle_state: Initialization Complete 2023-03-04 01:43:26 [INFO] idaes.init.fs.S101: Initialization Step 2 Complete: optimal - Optimal Solution Found 2023-03-04 01:43:26 [INFO] idaes.init.fs.C101.control_volume: Initialization Complete 2023-03-04 01:43:26 [INFO] idaes.init.fs.C101: Initialization Complete: optimal - Optimal Solution Found 2023-03-04 01:43:26 [INFO] idaes.init.fs.M101.mixed_state: Initialization Complete 2023-03-04 01:43:26 [INFO] idaes.init.fs.M101: Initialization Complete: optimal - Optimal Solution Found 2023-03-04 01:43:26 [INFO] idaes.init.fs.H101.control_volume: Initialization Complete 2023-03-04 01:43:26 [INFO] idaes.init.fs.H101: Initialization Complete: optimal - Optimal Solution Found 2023-03-04 01:43:26 [INFO] idaes.init.fs.R101.control_volume: Initialization Complete 2023-03-04 01:43:26 [INFO] idaes.init.fs.R101: Initialization Complete: optimal - Optimal Solution Found 2023-03-04 01:43:26 [INFO] idaes.init.fs.F101.control_volume: Initialization Complete 2023-03-04 01:43:26 [INFO] idaes.init.fs.F101: Initialization Complete: optimal - Optimal Solution Found 2023-03-04 01:43:26 [INFO] idaes.init.fs.S101.purge_state: Initialization Complete 2023-03-04 01:43:26 [INFO] idaes.init.fs.S101.recycle_state: Initialization Complete 2023-03-04 01:43:26 [INFO] idaes.init.fs.S101: Initialization Step 2 Complete: optimal - Optimal Solution Found 2023-03-04 01:43:26 [INFO] idaes.init.fs.C101.control_volume: Initialization Complete 2023-03-04 01:43:26 [INFO] idaes.init.fs.C101: Initialization Complete: optimal - Optimal Solution Found 2023-03-04 01:43:26 [INFO] idaes.init.fs.M101.mixed_state: Initialization Complete 2023-03-04 01:43:26 [INFO] idaes.init.fs.M101: Initialization Complete: optimal - Optimal Solution Found 2023-03-04 01:43:26 [INFO] idaes.init.fs.H101.control_volume: Initialization Complete 2023-03-04 01:43:26 [INFO] idaes.init.fs.H101: Initialization Complete: optimal - Optimal Solution Found 2023-03-04 01:43:26 [INFO] idaes.init.fs.R101.control_volume: Initialization Complete 2023-03-04 01:43:26 [INFO] idaes.init.fs.R101: Initialization Complete: optimal - Optimal Solution Found 2023-03-04 01:43:27 [INFO] idaes.init.fs.F101.control_volume: Initialization Complete 2023-03-04 01:43:27 [INFO] idaes.init.fs.F101: Initialization Complete: optimal - Optimal Solution Found 2023-03-04 01:43:27 [INFO] idaes.init.fs.S101.purge_state: Initialization Complete 2023-03-04 01:43:27 [INFO] idaes.init.fs.S101.recycle_state: Initialization Complete 2023-03-04 01:43:27 [INFO] idaes.init.fs.S101: Initialization Step 2 Complete: optimal - Optimal Solution Found 2023-03-04 01:43:27 [INFO] idaes.init.fs.C101.control_volume: Initialization Complete 2023-03-04 01:43:27 [INFO] idaes.init.fs.C101: Initialization Complete: optimal - Optimal Solution Found 2023-03-04 01:43:27 [INFO] idaes.init.fs.M101.mixed_state: Initialization Complete 2023-03-04 01:43:27 [INFO] idaes.init.fs.M101: Initialization Complete: optimal - Optimal Solution Found WARNING: Wegstein failed to converge in 3 iterations 2023-03-04 01:43:27 [INFO] idaes.init.fs.F102.control_volume: Initialization Complete 2023-03-04 01:43:27 [INFO] idaes.init.fs.F102: Initialization Complete: optimal - Optimal Solution Found Solving flowsheet... Ipopt 3.13.2: tol=1e-06 max_iter=5000 ****************************************************************************** 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...: 1031 Number of nonzeros in inequality constraint Jacobian.: 0 Number of nonzeros in Lagrangian Hessian.............: 934 Total number of variables............................: 340 variables with only lower bounds: 0 variables with lower and upper bounds: 146 variables with only upper bounds: 0 Total number of equality constraints.................: 340 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 6.60e+04 0.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0 1 0.0000000e+00 8.69e+03 1.42e+03 -1.0 2.00e+04 - 9.71e-01 4.67e-01H 1 2 0.0000000e+00 3.05e+03 1.56e+03 -1.0 1.60e+04 - 9.79e-01 4.90e-01h 1 3 0.0000000e+00 1.58e+03 1.55e+05 -1.0 1.41e+04 - 9.90e-01 4.99e-01h 1 4 0.0000000e+00 5.49e+02 8.87e+08 -1.0 8.43e+03 - 1.00e+00 9.57e-01h 1 5 0.0000000e+00 4.25e+03 2.87e+10 -1.0 8.02e+02 - 1.00e+00 9.90e-01h 1 6 0.0000000e+00 2.25e+03 1.51e+10 -1.0 8.39e+00 - 1.00e+00 1.00e+00h 1 7 0.0000000e+00 2.27e+01 1.40e+08 -1.0 2.45e-03 - 1.00e+00 1.00e+00f 1 8 0.0000000e+00 2.45e-03 1.23e+04 -1.0 2.38e-05 - 1.00e+00 1.00e+00h 1 9 0.0000000e+00 2.24e-08 3.40e-01 -2.5 1.06e-07 - 1.00e+00 1.00e+00h 1 Cannot recompute multipliers for feasibility problem. Error in eq_mult_calculator Number of Iterations....: 9 (scaled) (unscaled) Objective...............: 0.0000000000000000e+00 0.0000000000000000e+00 Dual infeasibility......: 2.8284425679755700e+05 2.8284425679755700e+05 Constraint violation....: 5.8207660913467407e-11 2.2351741790771484e-08 Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00 Overall NLP error.......: 5.8207660913467407e-11 2.8284425679755700e+05 Number of objective function evaluations = 11 Number of objective gradient evaluations = 10 Number of equality constraint evaluations = 11 Number of inequality constraint evaluations = 0 Number of equality constraint Jacobian evaluations = 10 Number of inequality constraint Jacobian evaluations = 0 Number of Lagrangian Hessian evaluations = 9 Total CPU secs in IPOPT (w/o function evaluations) = 0.014 Total CPU secs in NLP function evaluations = 0.001 EXIT: Optimal Solution Found. Complete.
IDAES Process Costing Framework¶
IDAES provides a library of capital costing correlations based on those in the following source:
Process and Product Design Principles: Synthesis, Analysis, and Evaluation. Seider, Seader, Lewin, Windagdo, 3rd Ed. John Wiley and Sons Chapter 22. Cost Accounting and Capital Cost Estimation 22.2 Cost Indexes and Capital Investment.
Currently, IDAES supports calculation of capital costing for a wide array of unit operations, vesseel sizing and material properties, and specific unit options such as column tray types and heat exchanger configurations. Users may find further information on specific costing methods and options in the IDAES Process Costing Framework documentation (link pending).
Add Operating Cost Equations¶
Before adding capital costing blocks, we will add operating cost equations taken from the basic HDA with Flash and the HDA with Distillation examples. The examples assume constant cooling and heating coefficients over an annual cost basis. The IDAES Generic Costing Framework does not currently support variable cost calculations.
# Required imports
from pyomo.environ import Expression
# Operating costs for HDA with second flash (model m)
m.fs.cooling_cost = Expression(expr=0.212e-7 * (-m.fs.F101.heat_duty[0]) +
0.212e-7 * (-m.fs.R101.heat_duty[0]))
m.fs.heating_cost = Expression(expr=2.2e-7 * m.fs.H101.heat_duty[0] +
1.9e-7 * m.fs.F102.heat_duty[0])
m.fs.operating_cost = Expression(expr=(3600 * 24 * 365 *
(m.fs.heating_cost +
m.fs.cooling_cost)))
Add Capital Costing¶
Below, we will add add capital costing blocks to the imported flowsheets and evaluate the economic impact of replacing the second Flash with a Distillation column. First, let's import and define the main flowsheet costing block:
# Import costing methods - classes, heaters, vessels, compressors, columns
from idaes.models.costing.SSLW import (
SSLWCosting,
SSLWCostingData,
)
from idaes.core import UnitModelCostingBlock
# Costing block
m.fs.costing = SSLWCosting()
Next, we will build the relevant costing blocks for the equipment we wish to cost. Note how the costing block, methods and flags are passed as arguments in the costing block call itself. Each unit model will have a single costing block, but each flowsheet model (m and n) will also have a single costing block for flowsheet-level properties.
Users should note that IDAES costing methods support a wide array of heating sources (e.g. fired, steam boiler, hot water) and do not support direct capital costing of coolers. If users wish to cost Heater units acting as coolers, it is necessary to cost a "dummy" 0D shell and tube exchanger with appropriate aliased hot stream properties and proper cooling water properties. This is not demonstrated here, as the HDA examples take advantage of Flash and Condenser operations to recover liquid product.
Capital costing is independent of unit model connections, and building cost equations may be done piecewise in this fashion. Default options are passed explicitly to demonstrate proper syntax and usage. Now that all required properties are defined, let's cost our models connecting costing blocks, methods and unit models in each flowsheet.
Flexibility of Costing Block Definitions¶
IDAES supports many ways to define batches of costing blocks, and several are shown in the example. Users may employ whichever method fits their modeling needs for explicit or concise code. In the code below, note how the unit model itself is never passed to the costing method; when the full model is executed, the costing block will automatically connect its parent block with child equation blocks.
Compressor
unit models with isothermal or adiabatic thermodynamics are too simple to cost and are therefore excluded from the economic analysis.
Let's define costing for the heater unit:
from idaes.models.costing.SSLW import (
HeaterMaterial,
HeaterSource,
)
# Costing for heater - m.fs.H101
m.fs.H101.costing = UnitModelCostingBlock(
flowsheet_costing_block=m.fs.costing,
costing_method=SSLWCostingData.cost_fired_heater,
costing_method_arguments={
"material_type": HeaterMaterial.CarbonSteel,
"heat_source": HeaterSource.Fuel,
}
)
The costing module provides a unit_mapping
dictionary linking generic unit model classes with recommended costing methods. In this example, StoichiometricReactor and Flash vessels utilize different vessel costing methods with similar arguments. The diameter and length attributes need to exist in order to cost vessel sizing and material requirements, and we add them if they don't exist already. The unit_mapping
method provides an opportunity to automatically select the correct vessel orientation (vertical or horizontal) based on the unit type; passing a StoichiometricReactor
or PFR
class object will call the cost_horizontal_vessel
method, while passing a Flash
or CSTR
class object will call the cost_vertical_vessel
method.
All vessels are assigned costing succintly via a loop below - users may define each block individually if desired:
from idaes.models.costing.SSLW import (
VesselMaterial,
TrayType,
TrayMaterial,
)
from idaes.core.util.constants import Constants
from pyomo.environ import Var, Constraint, units as pyunits, Param, value
from idaes.models.unit_models import StoichiometricReactor, Flash
# Map unit models to unit classes
# Will pass to unit_mapping which calls costing methods based on unit class
unit_class_mapping = {m.fs.R101: StoichiometricReactor,
m.fs.F101: Flash,
m.fs.F102: Flash,}
# Costing for vessels - m.fs.R101, m.fs.F101, m.fs.F102
# Loop over units
for unit in [m.fs.R101, m.fs.F101, m.fs.F102]:
# Get correct unit class for unit model
unit_class = unit_class_mapping[unit]
# Add dimension variables and constraint if they don't exist
if not hasattr(unit, "diameter"):
unit.diameter = Var(initialize=1, units=pyunits.m)
if not hasattr(unit, "length"):
unit.length = Var(initialize=1, units=pyunits.m)
if hasattr(unit, "volume"): # if volume exists, set diameter from volume
unit.volume_eq = Constraint(expr=unit.volume[0] == unit.length * unit.diameter**2 * 0.25 * Constants.pi)
else: # fix diameter directly
unit.diameter.fix(0.2214 * pyunits.m)
# Either way, fix L/D to calculate L from D
unit.L_over_D = Constraint(expr=unit.length == 3 * unit.diameter)
# Define vessel costing
unit.costing = UnitModelCostingBlock(
flowsheet_costing_block=unit.parent_block().costing, # e.g. m.fs.R101.costing
costing_method=SSLWCostingData.unit_mapping[unit_class], # e.g. cost_vertical_vessel()
costing_method_arguments={
"material_type": VesselMaterial.CarbonSteel,
"shell_thickness": 1.25 * pyunits.inch,
}
)
Solve Flowsheet Costing Blocks¶
Now, we may solve the full flowsheet for all costing blocks:
# Define solver
from idaes.core.solvers import get_solver
solver = get_solver()
# Check that the degrees of freedom are zero
from idaes.core.util.model_statistics import degrees_of_freedom
assert degrees_of_freedom(m) == 0
# Check physical units consistency, solve and check solver status
from pyomo.environ import TerminationCondition
from pyomo.util.check_units import assert_units_consistent
assert_units_consistent(m)
results = solver.solve(m, tee=True, symbolic_solver_labels=True)
assert results.solver.termination_condition == TerminationCondition.optimal
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...: 1067 Number of nonzeros in inequality constraint Jacobian.: 0 Number of nonzeros in Lagrangian Hessian.............: 942 Total number of variables............................: 358 variables with only lower bounds: 15 variables with lower and upper bounds: 146 variables with only upper bounds: 0 Total number of equality constraints.................: 358 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 8.81e+04 5.49e+03 -1.0 9.20e+04 - 2.62e-01 1.20e-01h 1 2 0.0000000e+00 6.34e+04 2.46e+03 -1.0 8.10e+04 - 5.83e-01 2.87e-01h 1 3 0.0000000e+00 3.40e+04 5.65e+03 -1.0 5.80e+04 - 5.72e-01 4.86e-01h 1 4 0.0000000e+00 2.33e+04 5.28e+05 -1.0 3.01e+04 - 7.04e-01 5.09e-01h 1 5 0.0000000e+00 1.13e+04 2.99e+09 -1.0 1.49e+04 - 8.02e-01 9.65e-01h 1 6 0.0000000e+00 5.50e+03 1.54e+09 -1.0 8.18e+02 - 9.90e-01 6.10e-01h 1 7 0.0000000e+00 5.23e+03 1.44e+09 -1.0 3.84e+02 - 9.92e-01 5.07e-02h 1 8 0.0000000e+00 5.22e+03 1.44e+09 -1.0 3.70e+02 - 1.00e+00 5.85e-04h 1 9r 0.0000000e+00 5.22e+03 1.00e+03 2.3 0.00e+00 - 0.00e+00 3.66e-07R 5 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 10r 0.0000000e+00 6.47e+04 1.34e+04 2.3 1.78e+05 - 2.30e-03 1.09e-03f 1 11r 0.0000000e+00 3.63e+04 4.53e+04 2.3 7.67e+03 - 7.32e-02 3.13e-03f 1 12r 0.0000000e+00 3.93e+04 4.20e+04 2.3 6.14e+03 - 6.79e-02 1.14e-01f 1 13r 0.0000000e+00 3.99e+04 2.20e+04 2.3 6.23e+02 - 4.86e-01 3.05e-01f 1 14r 0.0000000e+00 5.25e+04 1.67e+04 2.3 1.89e+03 - 2.37e-01 1.00e+00f 1 15r 0.0000000e+00 5.74e+04 8.19e+02 2.3 1.50e+03 - 9.73e-01 7.07e-01f 1 16r 0.0000000e+00 5.79e+04 7.74e+02 0.9 5.93e+02 - 9.40e-01 5.52e-01f 1 17r 0.0000000e+00 5.77e+04 4.91e+03 0.9 4.87e+02 - 9.04e-01 4.24e-01f 1 18r 0.0000000e+00 6.06e+04 1.62e+03 0.9 5.45e+02 - 1.00e+00 6.92e-01f 1 19r 0.0000000e+00 6.61e+04 2.32e+02 0.9 2.04e+02 - 1.00e+00 9.52e-01f 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 20r 0.0000000e+00 7.62e+04 1.03e+03 0.2 1.71e+02 - 8.32e-01 8.32e-01f 1 21r 0.0000000e+00 7.57e+04 8.29e+01 0.2 6.76e+02 - 1.00e+00 1.00e+00f 1 22r 0.0000000e+00 7.23e+04 4.74e+01 -0.5 1.37e+02 - 1.00e+00 9.01e-01f 1 23r 0.0000000e+00 5.71e+04 3.30e+01 -0.5 7.91e+02 - 1.00e+00 1.00e+00f 1 24r 0.0000000e+00 5.85e+04 6.01e-01 -0.5 1.26e+02 - 1.00e+00 1.00e+00h 1 25r 0.0000000e+00 4.88e+04 1.78e+02 -2.8 4.53e+02 - 9.80e-01 7.63e-01f 1 26r 0.0000000e+00 2.12e+04 2.63e+02 -2.8 2.29e+03 - 4.88e-01 2.16e-01f 1 27r 0.0000000e+00 2.05e+04 5.27e+02 -2.8 3.17e+03 - 1.00e+00 3.36e-01f 1 28r 0.0000000e+00 2.00e+04 7.36e+02 -2.8 3.93e+03 - 3.18e-02 1.06e-01f 1 29r 0.0000000e+00 7.19e+03 7.22e+02 -2.8 4.19e+03 - 1.00e+00 3.22e-01f 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 30r 0.0000000e+00 8.95e+03 2.50e+03 -2.8 6.03e+02 - 1.00e+00 7.72e-01f 1 31r 0.0000000e+00 5.53e+03 1.10e+03 -2.8 1.82e+03 - 7.92e-01 3.40e-01f 1 32r 0.0000000e+00 3.65e+03 7.38e+02 -2.8 1.83e+03 - 6.45e-01 5.76e-01f 1 33r 0.0000000e+00 4.19e+02 7.29e+02 -2.8 8.86e+02 - 6.85e-01 8.10e-01f 1 34r 0.0000000e+00 2.59e+02 3.88e+02 -2.8 8.93e+01 - 1.00e+00 1.00e+00f 1 35r 0.0000000e+00 2.37e+02 3.62e+02 -2.8 8.07e+00 - 1.00e+00 1.00e+00h 1 36r 0.0000000e+00 2.37e+02 1.60e+00 -2.8 2.49e+00 - 1.00e+00 1.00e+00h 1 37r 0.0000000e+00 2.37e+02 3.54e-02 -2.8 6.73e-02 - 1.00e+00 1.00e+00h 1 38r 0.0000000e+00 2.40e+02 6.94e+02 -4.3 1.16e+02 - 1.00e+00 5.13e-01f 1 39r 0.0000000e+00 1.40e+02 1.92e+03 -4.3 6.99e+01 - 9.98e-01 4.06e-01f 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 40r 0.0000000e+00 8.34e+01 2.11e+03 -4.3 4.21e+01 - 1.00e+00 6.67e-01f 1 41r 0.0000000e+00 2.54e+01 6.31e+02 -4.3 2.23e+01 - 7.35e-01 1.00e+00f 1 42 0.0000000e+00 2.26e+01 1.16e+07 -1.0 4.76e+01 - 9.90e-01 1.25e-01h 4 43 0.0000000e+00 2.00e+01 5.45e+06 -1.0 4.09e+01 - 9.90e-01 1.25e-01h 4 44 0.0000000e+00 1.59e+01 4.92e+06 -1.0 3.49e+01 - 9.91e-01 2.50e-01h 3 45 0.0000000e+00 1.23e+01 4.33e+06 -1.0 2.41e+01 - 1.00e+00 2.50e-01h 3 46 0.0000000e+00 7.56e+00 5.75e+06 -1.0 1.54e+01 - 1.00e+00 5.00e-01h 2 47 0.0000000e+00 3.06e+01 1.36e+07 -1.0 1.62e+00 - 1.00e+00 1.00e+00h 1 48 0.0000000e+00 3.18e-04 1.31e+03 -1.0 3.34e+01 - 1.00e+00 1.00e+00h 1 49 0.0000000e+00 2.98e-08 2.61e-01 -3.8 9.06e-08 - 1.00e+00 1.00e+00h 1 Cannot recompute multipliers for feasibility problem. Error in eq_mult_calculator Number of Iterations....: 49 (scaled) (unscaled) Objective...............: 0.0000000000000000e+00 0.0000000000000000e+00 Dual infeasibility......: 1.5041500809430196e+04 1.5041500809430196e+04 Constraint violation....: 2.9103830456733704e-11 2.9802322387695312e-08 Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00 Overall NLP error.......: 2.9103830456733704e-11 1.5041500809430196e+04 Number of objective function evaluations = 74 Number of objective gradient evaluations = 19 Number of equality constraint evaluations = 74 Number of inequality constraint evaluations = 0 Number of equality constraint Jacobian evaluations = 51 Number of inequality constraint Jacobian evaluations = 0 Number of Lagrangian Hessian evaluations = 49 Total CPU secs in IPOPT (w/o function evaluations) = 0.066 Total CPU secs in NLP function evaluations = 0.006 EXIT: Optimal Solution Found.
For comparison, we will call and build the HDA flowsheet replacing the second Flash
with a TrayColumn
distillation unit model. The flowsheet costing occurs in the external script hda_flowsheets_for_costing_notebook.py
and is not shown here:
from pyomo.common.log import LoggingIntercept
import logging
from io import StringIO
stream = StringIO()
with LoggingIntercept(stream, "idaes", logging.WARNING):
# Source file for prebuilt flowsheets
from hda_flowsheets_for_costing_notebook import hda_with_distillation
# Build hda model with distillation column and return model object
n = hda_with_distillation(tee=False)
Building flowsheet... Setting inputs... Initializing flowsheet... Limiting Wegstein tear to 3 iterations to obtain initial solution, if not converged IPOPT will pick up and continue. WARNING: Wegstein failed to converge in 3 iterations Solving flowsheet... WARNING: model contains export suffix 'fs.H102.control_volume.properties_out[0.0].scaling_factor' that contains 1 component keys that are not exported as part of the NL file. Skipping. WARNING: model contains export suffix 'fs.H102.control_volume.scaling_factor' that contains 1 component keys that are not exported as part of the NL file. Skipping. WARNING: model contains export suffix 'fs.F101.control_volume.properties_out[0.0].scaling_factor' that contains 26 component keys that are not exported as part of the NL file. Skipping. WARNING: model contains export suffix 'fs.F101.control_volume.properties_in[0.0].scaling_factor' that contains 25 component keys that are not exported as part of the NL file. Skipping. WARNING: model contains export suffix 'fs.F101.control_volume.scaling_factor' that contains 1 component keys that are not exported as part of the NL file. Skipping. WARNING: model contains export suffix 'fs.R101.control_volume.properties_out[0.0].scaling_factor' that contains 25 component keys that are not exported as part of the NL file. Skipping. WARNING: model contains export suffix 'fs.R101.control_volume.properties_in[0.0].scaling_factor' that contains 25 component keys that are not exported as part of the NL file. Skipping. WARNING: model contains export suffix 'fs.R101.control_volume.scaling_factor' that contains 2 component keys that are not exported as part of the NL file. Skipping. WARNING: model contains export suffix 'fs.H101.control_volume.properties_out[0.0].scaling_factor' that contains 26 component keys that are not exported as part of the NL file. Skipping. WARNING: model contains export suffix 'fs.H101.control_volume.properties_in[0.0].scaling_factor' that contains 25 component keys that are not exported as part of the NL file. Skipping. Adding distillation column and resolving flowsheet... WARNING: model contains export suffix 'fs.D101.condenser.control_volume.properties_out[0.0].scaling_factor' that contains 1 component keys that are not exported as part of the NL file. Skipping. WARNING: model contains export suffix 'fs.H102.control_volume.properties_out[0.0].scaling_factor' that contains 1 component keys that are not exported as part of the NL file. Skipping. WARNING: model contains export suffix 'fs.H102.control_volume.scaling_factor' that contains 1 component keys that are not exported as part of the NL file. Skipping. WARNING: model contains export suffix 'fs.F101.control_volume.properties_out[0.0].scaling_factor' that contains 26 component keys that are not exported as part of the NL file. Skipping. WARNING: model contains export suffix 'fs.F101.control_volume.properties_in[0.0].scaling_factor' that contains 25 component keys that are not exported as part of the NL file. Skipping. WARNING: model contains export suffix 'fs.F101.control_volume.scaling_factor' that contains 1 component keys that are not exported as part of the NL file. Skipping. WARNING: model contains export suffix 'fs.R101.control_volume.properties_out[0.0].scaling_factor' that contains 25 component keys that are not exported as part of the NL file. Skipping. WARNING: model contains export suffix 'fs.R101.control_volume.properties_in[0.0].scaling_factor' that contains 25 component keys that are not exported as part of the NL file. Skipping. WARNING: model contains export suffix 'fs.R101.control_volume.scaling_factor' that contains 2 component keys that are not exported as part of the NL file. Skipping. WARNING: model contains export suffix 'fs.H101.control_volume.properties_out[0.0].scaling_factor' that contains 26 component keys that are not exported as part of the NL file. Skipping. WARNING: model contains export suffix 'fs.H101.control_volume.properties_in[0.0].scaling_factor' that contains 25 component keys that are not exported as part of the NL file. Skipping. Complete. WARNING: model contains export suffix 'fs.D101.condenser.control_volume.properties_out[0.0].scaling_factor' that contains 1 component keys that are not exported as part of the NL file. Skipping. WARNING: model contains export suffix 'fs.H102.control_volume.properties_out[0.0].scaling_factor' that contains 1 component keys that are not exported as part of the NL file. Skipping. WARNING: model contains export suffix 'fs.H102.control_volume.scaling_factor' that contains 1 component keys that are not exported as part of the NL file. Skipping. WARNING: model contains export suffix 'fs.F101.control_volume.properties_out[0.0].scaling_factor' that contains 26 component keys that are not exported as part of the NL file. Skipping. WARNING: model contains export suffix 'fs.F101.control_volume.properties_in[0.0].scaling_factor' that contains 25 component keys that are not exported as part of the NL file. Skipping. WARNING: model contains export suffix 'fs.F101.control_volume.scaling_factor' that contains 1 component keys that are not exported as part of the NL file. Skipping. WARNING: model contains export suffix 'fs.R101.control_volume.properties_out[0.0].scaling_factor' that contains 25 component keys that are not exported as part of the NL file. Skipping. WARNING: model contains export suffix 'fs.R101.control_volume.properties_in[0.0].scaling_factor' that contains 25 component keys that are not exported as part of the NL file. Skipping. WARNING: model contains export suffix 'fs.R101.control_volume.scaling_factor' that contains 2 component keys that are not exported as part of the NL file. Skipping. WARNING: model contains export suffix 'fs.H101.control_volume.properties_out[0.0].scaling_factor' that contains 26 component keys that are not exported as part of the NL file. Skipping. WARNING: model contains export suffix 'fs.H101.control_volume.properties_in[0.0].scaling_factor' that contains 25 component keys that are not exported as part of the NL file. Skipping.
Results Comparison and Visualization¶
For the two flowsheets above, let's sum the total operating and capital costs of each scenario. We will display overall process economics results and compare the two flowsheets:
# Imports and data gathering
from matplotlib import pyplot as plt
plt.style.use('dark_background') # if using browser in dark mode, uncomment this line
import numpy as np
import pandas as pd
# Automatically get units that we costed - this will exclude C101 for both flowsheets
two_flash_unitlist = [getattr(m.fs, unit) for unit in dir(m.fs)
if hasattr(getattr(m.fs, unit), "costing")]
distillation_unitlist = [getattr(n.fs, unit) for unit in dir(n.fs)
if hasattr(getattr(n.fs, unit), "costing")]
# Compare equipment purchase costs (actual capital costs)
two_flash_capcost = {unit.name: value(unit.costing.capital_cost/1e3)
for unit in two_flash_unitlist}
distillation_capcost = {unit.name: value(unit.costing.capital_cost/1e3)
for unit in distillation_unitlist}
two_flash_capdf = pd.DataFrame(list(two_flash_capcost.items()),
columns=['Equipment','Two Flash']).set_index('Equipment')
distillation_capdf = pd.DataFrame(list(distillation_capcost.items()),
columns=['Equipment','Distillation']).set_index('Equipment')
# Add dataframes, merge same indices, replace NaNs with 0s, and transpose
capcosts = two_flash_capdf.add(distillation_capdf, fill_value=0).fillna(0).transpose()
# Sort according to an easier order to view
capcosts = capcosts[['fs.H101', 'fs.R101', 'fs.F101', 'fs.F102', 'fs.D101', 'fs.H102']]
print('Costs in $1000:')
display(capcosts) # view dataframe before plotting
capplot = capcosts.plot(kind='bar', stacked=True, title='HDA Total Capital Costs', ylabel='$1000')
Costs in $1000:
Equipment | fs.H101 | fs.R101 | fs.F101 | fs.F102 | fs.D101 | fs.H102 |
---|---|---|---|---|---|---|
Distillation | 13.283334 | 13.552294 | 7.781332 | 0.000000 | 7.60743 | 3.32944 |
Two Flash | 13.283334 | 10.755878 | 7.781332 | 7.781332 | 0.00000 | 0.00000 |
# Compare operating costs (per year)
two_flash_opcost = {'cooling': value(3600 * 24 * 365 * m.fs.cooling_cost/1e3),
'heating': value(3600 * 24 * 365 * m.fs.heating_cost/1e3)}
distillation_opcost = {'cooling': value(3600 * 24 * 365 * n.fs.cooling_cost/1e3),
'heating': value(3600 * 24 * 365 * n.fs.heating_cost/1e3)}
two_flash_opdf = pd.DataFrame(list(two_flash_opcost.items()),
columns=['Utilities','Two Flash']).set_index('Utilities')
distillation_opdf = pd.DataFrame(list(distillation_opcost.items()),
columns=['Utilities','Distillation']).set_index('Utilities')
# Add dataframes, merge same indices, replace NaNs with 0s, and transpose
opcosts = two_flash_opdf.add(distillation_opdf, fill_value=0).fillna(0).transpose()
print('Costs in $1000:')
display(opcosts) # view dataframe before plotting
opplot = opcosts.plot(kind='bar', stacked=True, title='HDA Operating Costs', ylabel='$1000/year')
Costs in $1000:
Utilities | cooling | heating |
---|---|---|
Distillation | 60.192008 | 367.404722 |
Two Flash | 47.028803 | 372.093536 |
# Compare total costs (capital costs and operating costs)
two_flash_totcost = {'capital': sum(two_flash_capcost[idx] for idx in two_flash_capcost),
'operating': value(m.fs.operating_cost)/1e3}
distillation_totcost = {'capital': sum(distillation_capcost[idx] for idx in distillation_capcost),
'operating': value(n.fs.operating_cost)/1e3}
two_flash_totdf = pd.DataFrame(list(two_flash_totcost.items()),
columns=['Costs','Two Flash']).set_index('Costs')
distillation_totdf = pd.DataFrame(list(distillation_totcost.items()),
columns=['Costs','Distillation']).set_index('Costs')
# Add dataframes, merge same indices, replace NaNs with 0s, and transpose
totcosts = two_flash_totdf.add(distillation_totdf, fill_value=0).fillna(0).transpose()
print('Costs in $1000:')
display(totcosts) # view dataframe before plotting
totplot = totcosts.plot(kind='bar', stacked=True, title='HDA Total Plant Cost (TPC)', ylabel='$1000/year')
Costs in $1000:
Costs | capital | operating |
---|---|---|
Distillation | 45.553829 | 427.596731 |
Two Flash | 39.601876 | 419.122339 |
Finally, let's compare the total costs on a production basis. This will account for the greater efficiency provided by the distillation column relative to the less-expensive second flash unit:
two_flash_cost = value(1e3*sum(two_flash_totcost[idx] for idx in two_flash_totcost))
two_flash_prod = value(m.fs.F102.vap_outlet.flow_mol_phase_comp[0, "Vap", "benzene"] * 365 * 24 * 3600)
distillation_cost = value(1e3*sum(distillation_totcost[idx] for idx in distillation_totcost))
distillation_prod = value(n.fs.D101.condenser.distillate.flow_mol[0] * 365 * 24 * 3600)
print(f"Two flash case over one year: ${two_flash_cost/1e3:0.0f}K / {two_flash_prod/1e3:0.0f} kmol benzene = ${two_flash_cost/(two_flash_prod/1e3):0.2f} per kmol benzene produced")
print(f"Distillation case over one year: ${distillation_cost/1e3:0.0f}K / {distillation_prod/1e3:0.0f} kmol benzene = ${distillation_cost/(distillation_prod/1e3):0.2f} per kmol benzene produced")
Two flash case over one year: $459K / 4477 kmol benzene = $102.45 per kmol benzene produced Distillation case over one year: $473K / 5108 kmol benzene = $92.63 per kmol benzene produced
Summary¶
In this example, IDAES Process Costing Framework methods were applied to two HDA flowsheets for capital cost estimation. The costing blocks calls showcased multiple methods to define unit costing, demonstrating the flexibility and best practice of the costing framework. In the basic examples, the two-flash HDA did not include costing and the distillation HDA estimated a reactor capital cost comprising 3.3% of the total plant cost (TPC). With more rigorous costing, IDAES obtained total capital costs of 8.5% TPC (two flash HDA) and 9.6% (distillation HDA) and better modeled the impact of equipment cost on process economics. As printed above, the IDAES Process Costing Framework confirmed that replacing the second flash drum with a distillation column results in increased equipment costs, increased production and decreased cost per unit product.