Show code cell content
###############################################################################
# The Institute for the Design of Advanced Energy Systems Integrated Platform
# Framework (IDAES IP) was produced under the DOE Institute for the
# Design of Advanced Energy Systems (IDAES).
#
# Copyright (c) 2018-2023 by the software owners: The Regents of the
# University of California, through Lawrence Berkeley National Laboratory,
# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon
# University, West Virginia University Research Corporation, et al.
# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md
# for full copyright and license information.
###############################################################################
HDA Flowsheet Costing#
Maintainer: Brandon Paul
Author: Brandon Paul
Updated: 2023-06-01
Note#
This example 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.pyidaes.models.properties.activity_coeff_models.BTX_activity_coeff_VLEhda_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...
WARNING: Params with units must be mutable. Converting Param
'fs.thermo_params.pressure_crit' to mutable.
WARNING: Params with units must be mutable. Converting Param
'fs.thermo_params.temperature_crit' to mutable.
WARNING: Params with units must be mutable. Converting Param
'fs.thermo_params.mw_comp' to mutable.
WARNING: Params with units must be mutable. Converting Param
'fs.thermo_params.dens_liq_param_1' to mutable.
WARNING: Params with units must be mutable. Converting Param
'fs.thermo_params.dens_liq_param_2' to mutable.
WARNING: Params with units must be mutable. Converting Param
'fs.thermo_params.dens_liq_param_3' to mutable.
WARNING: Params with units must be mutable. Converting Param
'fs.thermo_params.dens_liq_param_4' to mutable.
WARNING: Params with units must be mutable. Converting Param
'fs.thermo_params.temperature_boil' to mutable.
WARNING: Params with units must be mutable. Converting Param
'fs.thermo_params.cp_ig_1' to mutable.
WARNING: Params with units must be mutable. Converting Param
'fs.thermo_params.cp_ig_2' to mutable.
WARNING: Params with units must be mutable. Converting Param
'fs.thermo_params.cp_ig_3' to mutable.
WARNING: Params with units must be mutable. Converting Param
'fs.thermo_params.cp_ig_4' to mutable.
WARNING: Params with units must be mutable. Converting Param
'fs.thermo_params.cp_ig_5' to mutable.
WARNING: Params with units must be mutable. Converting Param
'fs.thermo_params.pressure_sat_coeff_A' to mutable.
WARNING: Params with units must be mutable. Converting Param
'fs.thermo_params.pressure_sat_coeff_B' to mutable.
WARNING: Params with units must be mutable. Converting Param
'fs.thermo_params.pressure_sat_coeff_C' to mutable.
WARNING: Params with units must be mutable. Converting Param
'fs.thermo_params.dh_vap' to mutable.
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
2024-05-09 18:33:54 [INFO] idaes.init.fs.H101.control_volume: Initialization Complete
---------------------------------------------------------------------------
ApplicationError Traceback (most recent call last)
Cell In[2], line 5
2 from hda_flowsheets_for_costing_notebook import hda_with_flash
4 # Build hda model with second flash unit and return model object
----> 5 m = hda_with_flash(tee=True)
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/checkouts/latest/idaes_examples/notebooks/docs/flowsheets/hda_flowsheets_for_costing_notebook.py:221, in hda_with_flash(tee)
218 def function(unit):
219 unit.initialize(outlvl=outlvl)
--> 221 seq.run(m, function)
223 # solve model
224 print('Solving flowsheet...')
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/network/decomposition.py:296, in SequentialDecomposition.run(self, model, function)
293 self.cache.clear()
295 try:
--> 296 return self._run_impl(model, function)
297 finally:
298 # Cleanup
299 self.cache.clear()
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/network/decomposition.py:317, in SequentialDecomposition._run_impl(self, model, function)
315 logger.info("Starting first pass run of network")
316 order = self.calculation_order(G)
--> 317 self.run_order(G, order, function, tset, use_guesses=True)
319 if not self.options["solve_tears"] or not len(tset):
320 # Not solving tears, we're done
321 end = time.time()
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/network/decomposition.py:407, in SequentialDecomposition.run_order(self, G, order, function, ignore, use_guesses)
404 self.load_guesses(guesses, port, fixed_ins)
405 self.load_values(port, default, fixed_ins, use_guesses)
--> 407 function(unit)
409 # free the inputs that were not already fixed
410 for var in fixed_ins:
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/checkouts/latest/idaes_examples/notebooks/docs/flowsheets/hda_flowsheets_for_costing_notebook.py:219, in hda_with_flash.<locals>.function(unit)
218 def function(unit):
--> 219 unit.initialize(outlvl=outlvl)
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/base/unit_model.py:540, in UnitModelBlockData.initialize(blk, *args, **kwargs)
537 c.deactivate()
539 # Remember to collect flags for fixed vars
--> 540 flags = blk.initialize_build(*args, **kwargs)
542 # If costing block exists, activate and initialize
543 for c in init_order:
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/base/unit_model.py:601, in UnitModelBlockData.initialize_build(blk, state_args, outlvl, solver, optarg)
598 # ---------------------------------------------------------------------
599 # Solve unit
600 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
--> 601 results = opt.solve(blk, tee=slc.tee)
603 init_log.info_high(
604 "Initialization Step 2 {}.".format(idaeslog.condition(results))
605 )
607 # ---------------------------------------------------------------------
608 # Release Inlet state
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/opt/base/solvers.py:534, in OptSolver.solve(self, *args, **kwds)
531 def solve(self, *args, **kwds):
532 """Solve the problem"""
--> 534 self.available(exception_flag=True)
535 #
536 # If the inputs are models, then validate that they have been
537 # constructed! Collect suffix names to try and import from solution.
538 #
539 from pyomo.core.base.block import BlockData
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/opt/solver/shellcmd.py:140, in SystemCallSolver.available(self, exception_flag)
138 if exception_flag:
139 msg = "No executable found for solver '%s'"
--> 140 raise ApplicationError(msg % self.name)
141 return False
142 return True
ApplicationError: No executable found for solver 'ipopt'
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 succinctly 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:
# Eefine solver
from idaes.core.solvers import get_solver
solver = get_solver()
# Check that the degrees of freedom is 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.24e-08 3.21e-01 -3.8 1.06e-07 - 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.5041240304748950e+04 1.5041240304748950e+04
Constraint violation....: 2.9103830456733704e-11 2.2351741790771484e-08
Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00
Overall NLP error.......: 2.9103830456733704e-11 1.5041240304748950e+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.061
Total CPU secs in NLP function evaluations = 0.003
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.